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ABSTRACT 


The  objective  in  preparing  this  report  is  to  provide  a  perspective 
on  the  signal  processing  problem  that  must  be  solved  in  order  to  produce 
a  useful  seismic  surveillance  system.  The  point-of-view  is  taken  that  the 
problem  is  nontrivial,  thereby  requiring  a  particularly  careful  development 
of  the  signal  processing  algorithms.  The  surveillance  system  is  faced  with 
the  problem  of  detecting  weak  signals  in  the  presence  of  incoherent  and 
coherent  noise  sources.  The  signals  that  are  sought  generally  have  an 
unknown  structure  except  for  some  important  qualitative  characteristics. 
There  qualitative  characteristics  are  identified  by  considering  the  propaga¬ 
tion  of  seismic  energy  in  terms  of  the  possible  wave  types,  propagation 
paths,  and  frequency  ranges  of  possible  signals. 

The  presentation  is  composed  of  three  major  sections.  In  Part  I, 
basic  physical  considerations  regarding  the  propagation  and  sensing  of 
seismic  energy  are  reviewed.  This  discussion  provides  physical  motiva¬ 
tion  for  the  development  of  suitable  signal  processing  algorithms.  The 
general  signal  processing  problem  is  defined  in  Part  II.  Array  process¬ 
ing  is  reviewed  and  fundamental  concepts  are  introduced  in  this  section. 
General  array  processing  algorithms  are  presented  in  Part  III.X  The 
presentation  in  this  section  is  intentionally  brief  as  the  algorithms  are 
motivated  and  developed  in  Parts  I  and  II. 

The  statement  of  the  algorithms  in  Part  III  is  very  general. 

It  is  the  contention  of  this  report  that  the  surveillance  problem  is  suf¬ 
ficiently  difficult  that  unnecessary  haste  must  be  avoided  in  developing 
specific  algorithms  of  limited  scope.  It  is  important  that  a  maximal 
amount  of  information  is  extracted  from  the  sensor  outputs  in  order  to 
achieve  desirable  performance  goals.  Thus,  the  general  structure  of 
the  algorithms  is  defined  to  serve  as  a  basis  for  the  evolution  of  spe¬ 
cific  algorithms  that  can  be  implemented  in  the  seismic  surveillance 
system.  By  proceeding  early  in  the  development  to  specific  algorithms 
based  upon  restrictive  assumptions,  one  can  be  misled  in  assessing  the 
potential  performance  of  the  system. 


ORINCON 


i 


INTRODUCTION  AND  OVERVIEW 


A  primary  objective  of  the  signal  processor  is  to  detect  seismic 
energy  from  an  unknown  source  and  to  estimate  the  location  of  the 
source.  Many  sources  are  typified  as  introducing  an  impulse-like  signal 
into  the  earth  (e.g.,  an  artillery  recoil  or  a  shell  blast).  Thus,  the 
input  signal  occurs  during  a  relatively  short  time  interval.  Because  of 
the  short  time  duration,  the  signal  will  contain  a  broad  range  of  fre¬ 
quencies.  Alternative  sources  of  energy  are  provided  by  moving  tar¬ 
gets  (e.g.,  tanks,  trucks,  aircraft).  These  sources  introduce  energy 
into  the  earth  over  longer  time  intervals  and  generally  can  be  expected 
to  define  signals  having  narrowband  characteristics.  Both  types  of 
sources  are  considered  in  the  following  discussion. 

The  earth  can  be  regarded  as  a  channel  that  transmits  the 
seismic  energy  generated  by  a  source  to  sensors  that  provide  infor¬ 
mation  for  the  signal  processor.  The  effect  of  the  earth  upon  the 
input  signal  is  very  complex.  A  single  source  can  generate  several 
signals  that  may  be  observed  by  the  sensors.  These  signals  repre¬ 
sent  the  different  types  of  waves  that  can  propagate  in  an  elastic 
medium  or  the  different  propagation  paths  that  are  possible  for  a 
specific  type  of  wave.  To  design  the  signal  processor,  the  character¬ 
istics  of  the  wave  motion  and  the  possible  propagation  paths  must  be 
understood  and  modeled. 

The  designer  of  the  signal  processor  must  be  concerned  with 
the  properties  of  the  signals  received  by  the  sensors.  Because  of  the 
complex  nature  of  the  transmission  channel,  the  signal  received  by 
the  sensors  is  substantially  different  from  the  input  signal  generated 
by  the  source.  To  reduce  the  effects  of  the  channel,  the  "sensor" 
typically  consists  of  an  array  of  three-axis  geophones  whose  outputs 
are  combined  in  an  appropriate  manner  by  the  signal  processor.  The 
algorithms  that  define  the  processor  must  consider  the  following  aspects. 


Signal  Bandwidth:  The  range  of  frequencies  which  the 
signal  is  expected  to  contain  must  be  established.  As  the 
bandwidth  increases,  the  time  duration  of  the  signal  diminishes. 

Signal  Power:  The  difficulty  of  the  detection  problem  increases 
as  the  signal-to-noise  ratio  decreases.  The  signal  power  dimin¬ 
ishes  as  the  distance  between  the  source  and  the  sensor  increases. 
Because  the  surveillance  system  is  intended  to  operate  at  rela¬ 
tively  long  range  (i.e.,  10  to  20  km),  the  signal  processor  must 
be  designed  to  operate  at  low  signal-to-noise  ratios  (SNR). 

Basic  Wave  Properties:  Different  types  of  waves  induce  different 
particle  motions  and  these  characteristics  may  be  detected  by  a 
three-axis  geophone.  Then,  the  type  of  wave  can  be  established. 
Because  different  waves  propagate  with  different  speeds,  this 
classifications  can  be  used  in  the  signal  processor  to  enhance 
detection  and  to  aid  in  source  localization. 

Signals  that  are  received  by  the  sensor  can  result  from 
waves  that  propagate  along  a  variety  of  paths  from  the  source 
to  the  sensor.  Characterization  of  the  propagation  path  is 
useful  for  several  reasons.  For  example,  a  multiple  reflection 
can  arrive  at  a  time  that  is  sufficiently  later  than  other  arrivals 
that  it  appears  to  be  a  new  signal.  Thus,  the  detector  may 
introduce  a  false  alarm  that  can  be  otherwise  avoided. 

Wavefront  Characteristics:  Arrays  of  geophones  must  be  used 
in  order  to  obtain  an  increase  in  the  effective  SNR  relative  to 
a  single  geophone.  But  the  geometrical  separation  of  the 
elements  of  an  array  can  introduce  the  requirement  for  character¬ 
izing  the  properties  of  the  v/avefront  that  is  sensed  by  the 
geophones.  Planar  or  spherical  wavefronts  may  describe  the 
propagation  of  many  waves.  However,  in  some  circumstances 
the  waves  can  exhibit  a  substantial  distortion  of  the  wavefront 
that  must  be  considered  by  the  signal  processor,  particularly 
for  inter-array  processing. 


Error  Sources;  The  analysis  of  wave  propagation  in  the  earth 
begins  with  several  idealized  assumptions.  These  assumptions 
lead  to  the  general  descriptions  of  signal  bandwidth,  signal 
power  characteristics,  wave  properties,  and  wavefront  behavior 
mentioned  above.  Deviations  from  the  basic  assumptions  will 
introduce  errors  in  tne  general  models  and  the  designer  of  the 
signal  processor  must  be  sensitive  to  effects  of  the  potential 
error  sources.  These  include  effects  such  as  dips  in  the  layer¬ 
ing  of  the  earth,  the  influence  of  the  weathered  layer  at  the 
surface  of  the  earth,  elevation  differences  between  the  source 
and  the  sensor,  and  the  irregularities  in  the  subsurface  compo¬ 
sition. 

The  signal  processing  concerns  that  have  been  mentioned  are 
addressed  by  defining  idealized  properties  for  the  earth,  developing 
mathematical  models  from  these  assumptions,  and  analyzing  the  conse¬ 
quences  of  the  models.  In  Part  I  the  assumptions  and  their  conse¬ 
quences  are  reviewed  to  obtain  the  basic  features  that  must  be  con¬ 
sidered  in  the  design  of  the  signal  processing  algorithms.  In  Section 
1.2,  the  basic  assumptions  are  stated  and  the  models  are  reviewed. 

The  types  of  waves  and  their  general  characteristics  are  discussed  in 
Section  I.  3.  The  potential  propagation  paths  are  described  in  Section 
1.4.  Some  basic  errors,  sources  and  their  effects  are  investigated  in 
Section  1.5. 

Before  beginning  the  discussion  of  the  properties  of  the  trans¬ 
mission  channel,  a  general  description  of  the  sensors  is  provided  in 
Section  1.1.  This  discussion  is  intended  to  be  qualitative  and  to 
establish  a  motivation  for  the  subsequent  description  of  the  signal 
characteristics  presented  in  Sections  1.2  through  1.5.  The  comments 
of  Sections  1.1  through  1.5  are  summarized  in  Section  1.6.  In  this 
concluding  section,  some  basic  modeling  and  signal  processing  conclusions 
are  stated  and  discussed. 


The  analysis  of  the  propagation  of  seismic  energy  leads  to  the 
identification  of  several  wave  types  and  a  variety  of  propagation  paths. 
Each  wave  induces  a  distinctive  particle  motion  that  must  be  measured 
by  the  sensors  of  the  seismic  surveillance  system.  Each  wave  travels 
with  some  average  velocity  and  may  be  attenuated  and  dispersed  along 
its  path.  These  and  other  considerations  lead  to  the  conclusion  that 
the  signal  received  by  a  geophone  from  a  particular  seismic  event  will 
be  a  composite  of  different  wave  types  that  have  traveled  along  several 
propagation  paths.  Consequently,  the  signal  can  be  expected  to  have 
a  very  complicated  and,  largely,  unpredictable  character. 

While  the  form  of  a  signal  must  be  regarded  as  unpredictable, 
some  general  characteristics  can  be  assumed.  For  an  impulsive  source 
(e.g.,  an  artillery  recoil  or  a  shell  blast),  the  wave  duration  will  be 
relatively  brief,  lasting  for  at  most  a  few  seconds  so  that  a  broad  band 
signal  is  received.  Geophysical  considerations  imply  that  detectable 
power  is  restricted  to  frequencies  ranging  from  a  few  hertz  to  100  hertz. 
The  upper  limit  of  the  frequency  band  generally  diminishes  as  the  dis¬ 
tance  from  the  source  to  the  sensor  increases.  The  short-time  duration 
of  the  signals  induced  by  impulse-like  sources  suggests  the  need  for 
time-domain  rather  than  frequency-domain  processing. 

Other  target  sources  can  induce  narrowband  signals.  For  example, 
truck  or  tank  motion  imparts  a  persistent  input  to  the  earth  that  one 
would  expect  to  have  a  narrowband  characterization.  For  narrowband 
signals,  frequency-domain  processing  is  appropriate.  These  observa¬ 
tions  indicate  that  care  must  be  taken  in  the  signal  processing  to  recog¬ 
nize  the  type  of  signal  that  is  sought  and  to  apply  appropriate  processing 
algorithms.  In  subsequent  discussion,  both  wideband  and  narrowband 
processing  algorithms  shall  be  discussed. 

The  basic  function  of  the  surveillance  system  is  to  detect  a 
signal,  whether  it  is  narrowband  or  wideband.  Because  the  signal-to- 
noise  ratio  (i.e.,  SNR)  can  be  expected  to  be  small,  the  detection 
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problem  is  difficult.  A  second  important  function  of  the  system  is  to 
determine  the  direction  from  which  the  detected  s.gnal  emanates.  The 
detection  and  direction-finding  requirements  imply  that  an  array  of 
seisors  should  be  deployed.  Because  of  the  time-varying  character  of 
the  noise,  adaptive  array  processing  methods  should  be  employed  to 
achieve  the  system  objective.'..  There  may  be  several  sources  of  coher¬ 
ent  energy  during  any  detection  interval.  Therefore,  the  directional 
noise  cancelling  capabilities  of  an  adaptive  array  must  be  employed  and 
should  be  particularly  useful.  Further,  the  data  from  an  array  can  be 
used  to  estimate  the  average  velocity  of  a  detected  signal  as  it  passes 
the  array  location.  The  velocitv  estimates  can  be  used,  for  example,  as 
one  feature  for  establishing  the  type  of  wave.  The  surveillance  system, 
also,  should  provide  information  regarding  the  location  of  the  source  of 
a  detected  signal.  This  requirement  implies  the  utilization  of  multiple 
arrays.  By  processing  the  outputs  from  more  than  one  array,  localiza¬ 
tion  and  tracking  information  can  be  deduced. 

Mathematical  models  that  are  useful  for  developing  signal  process¬ 
ing  algorithms  are  introduced  in  Part  II.  Basic  models  and  the  general 
problems  are  defined  in  Section  II.  1.  Some  fundamental  aspects  of  array 
processing  are  discussed  in  Section  II.  2.  In  particular,  the  array  pattern 
is  defined  and  discussed  as  it  provides  a  basic  tool  for  evaluating  the 
performance  of  a  particular  array  configuration.  The  array  processing 
problem  is  described  in  Section  II.  3  and  the  concept  of  array  directivity 
is  introduced.  The  conditions  for  arrays  having  maximal  directivity  are 
presented  and  discussed.  Then,  a  statistical  formulation  of  the  optimal 
array  processing  problem  is  defined  and  a  general  solution  is  described 
in  Section  II.  4.  Additional  aspects  of  array  processing  are  introduced 
in  this  section  also.  Adaptive  arrays  are  discussed;  optimal  detectors 
are  introduced;  and  narrowband  and  broadband  signals  are  considered. 

The  salient  features  of  the  seismic  signal  processing  problem, 
discussed  in  Parts  I  and  II,  arc  summarized  in  Part  III.  General  algo¬ 
rithms  that  are  suitable  for  establishing  the  feasibility  of  the  system 


from  experimental  and/or  simulated  data  are  presented  in  Part  III. 
One  can  proceed  directly  from  this  Introduction  to  Part  III  to  avoid 
the  lengthier  discussions  ot  the  physical  character  of  the  problem 
and  the  major  signal  processing  concerns.  Emphasis  in  Part  III  is 
placed  on  the  presentation  of  the  general  computations  that  must  be 
accomplished  to  obtain  the  best  performance  of  the  signal  processor. 
In  addition,  recommendations  are  provided  for  directions  that  should 
be  taken  for  the  detailed  development  of  the  signal  processor. 


PART  I 


BASIC  PHYSICAL  CONSIDERATIONS 


ABSTRACT 

Before  defining  algorithms  that  can  be  used  in  a  seismic  surveil 
lance  system  for  detection,  classification,  localization,  and  tracking  of 
relevant  objects  and  events,  it  is  necessary  to  examine  the  character 
of  the  signals  received  by  the  sensors  in  order  to  assess  the  signal 
processing  requirements.  The  fundamental  properties  concerning  the 
propagation  of  seismic  energy  are  reviewed  in  this  section.  The  dis¬ 
cussion  shall  focus  upon  those  aspects  that  are  most  germane  to  the 
development  of  effective  signal  processing  algorithms.  More  detailed 
expositions  are  provided,  for  example,  in  the  books  by  Dobrin  [1], 
Grant  and  West  ( 2\  ,  Parasnis  [31,  Telford,  et.  al.  [4],  and  Aki  and 
Richards  [5|.  This  presentation  draws  heavily  from  these  references. 
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1 .  1  General  Description  of  the  Sensors 

The  seismic  energy  transmitted  within  the  earth  by  a  specific 
source  is  detected  by  a  geophone.  The  output  of  this  sensor  gener¬ 
ally  can  be  regarded  as  being  proportional  to  the  velocity  of  the  motion 
of  the  earth  at  the  geophone  in  the  direction  of  its  sensitive  axis. 

The  natural  frequency  of  a  geophone  is  typically  less  than 
10  Hz  and  may  be  as  low  as  1  or  2  Hz.  The  instruments  exhibit 
a  flat  amplitude  response  and  linear  phase  for  frequencies  between 
the  natural  frequency  and,  possibly,  several  hundred  Hertz.  The 
range  of  possible  frequencies  in  each  type  of  seismic  wave  is  dis¬ 
cussed  below.  It  shall  be  assumed  that  the  geophones  are  sensitive 
to  the  range  of  potential  frequencies  that  may  be  encountered. 

Motions  in  both  horizontal  and  vertical  directions  are  possible 
and  can  provide  useful  information.  Consequently,  three  geophones 
can  be  mounted  together  to  form  a  three-axis  geophone.  Typically, 
the  sensitive  axes  are  aligned  to  be  orthogonal  to  one  another.  Then, 
the  three-axis  geophone  is  placed  on  or  near  the  surface  of  the 
earth  with  one  geophone  oriented  vertically  and  with  the  two 
remaining  geophones  located  in  the  horizontal  plane.  In  the  subse¬ 
quent  discussion,  reference  will  be  made  to  particle  motions  in  the 
horizontal  (i.e.  ,  X  and  Y  directions)  and  vertical  (i.e..  Z  direction) 
directions  in  implicit  reference  to  this  assumed  orientation. 

The  signal  received  by  a  geophone  generally  can  be  expected 
to  be  small  relative  to  the  noise.  To  enhance  the  signal,  arrays  of 
geophones  are  formed  whose  outputs,  when  combined  in  an  appro¬ 
priate  manner,  produce  a  greater  SNR  in  the  composite  signal  than 
in  the  elementary  signal  obtained  from  a  single  sensor.  The  array 
configuration  can  be  defined  in  a  variety  of  ways.  Questions  of  this 
nature  are  an  integral  part  of  the  design  of  the  signal  processor 
and  will  be  considered  later  in  the  discussion. 
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It  is  also  natural  to  consider  more  than  one  array.  In  other 
words,  it  is  desirable  to  consider  the  placement  of  arrays  at  widely 
separated  locations.  The  outputs  of  each  array  can  provide  a  sig¬ 
nificantly  different  view  of  a  source  to  assist  not  only  with  detection 
but  with  localization  and  tracking.  Because  of  the  existence  of 
arrays  of  sensors,  the  spatial  characteristics  of  the  propagation  of 
seismic  energy  must  be  considered.  However,  the  discussion  of  the 
signal  processing  problem  presented  in  this  report  is  restricted  to  a 
single  array.  Multi-array  processing,  while  an  important  topic,  is 
beyond  the  scope  of  the  current  effort. 
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Basic  Modeling  Assumptions 


The  mathematical  analysis  of  the  propagation  of  seismic  energy 
begins  by  considering  an  idealized  model  of  the  earth  that  can  be  used 
to  obtain  fund  •—  •  d  insights.  Then,  the  effect  of  deviations  from 
the  ideal  can  be  considered  to  obtain  an  understanding  of  the  propa¬ 
gation  for  more  realistic  models.  Fundamental  to  these  analyses  is 
the  assumption  that  the  earth  is  an  elastic  medium.  Consequently, 
one  is  concerned  primarily  with  the  nature  of  elastic-wave  propa¬ 
gation  .  This  assumption  imposes  a  linear  character  upon  the  channel 
characterization  that  is  important  from  the  point-of-view  of  the 
signal  processor. 

The  propagation  of  seismic  energy  occurs  through  particle 
deformations  that  travel  through  the  earth  at  velocities  that  are 
determined  by  the  elastic  properties  and  densities  of  material.  The 
deformations  are  analyzed  in  terms  of  stress  and  strain  (i.e.,  the 
forces  that  cause  the  deformations) .  The  analysis  leads  to  the  devel¬ 
opment  of  the  wave  equation  as  the  descriptor  of  the  particle  motion. 

The  elastic  nature  of  the  material  provides  the  conclusion  that 
there  are  two  basic  types  of  waves  that  must  be  considered, 
compressional  waves  and  shear  waves.  Generally  referred  to  as 
body  waves ,  the  propagation  of  these  waves  must  be  investigated 
to  establish  the  manner  in  which  the  seismic  energy  is  absorbed 
and  attenuated  by  the  material. 

Wave  propagation  is  studied  by  assuming  an  idealized  earth 
composed  of  horizontal  layers  of  materials  with  radically  different 
elastic  properties.  The  interface  between  layers  is  assumed  to  be 
described  by  a  discontinuous  change  in  the  elastic  properties. 

Within  a  layer,  it  is  commonly  assumed  that  the  material  is  homo¬ 
geneous  and  isotropic.  With  this  model,  the  propagation  paths  for 
the  body  waves  can  be  analyzed  in  terms  of  their  reflection  and 
refraction  characteristics  at  the  interface  boundaries. 


In  addition  to  the  compressional  and  shear  waves,  a  Rayleigh 
wave  travels  along  the  free  surface  of  the  earth.  As  discussed 
below,  the  Rayleigh  wave  is,  generally,  the  strongest  wave  gener¬ 
ated  by  a  compressional  source.  The  layering  of  the  earth  can 
permit  the  existence  of  another  surface  wave,  the  Love  wave.  The  Love 
wave  is  observed  when  a  low-velocity  layer  overlays  a  higher-speed 
layer.  Information  regarding  the  surface  or  weathered  layer  can  be 
obtained  from  the  Love  wave. 

Deviations  from  the  idealized  model  can  complicate  the 
analysis  in  a  significant  fashion.  Imperfectly  elastic,  non-isotropic, 
or  non-homogeneous  materials  must  be  considered.  Layering  generally 
will  not  be  perfectly  horizontal.  In  fact,  the  interface  may  not  be 
represented  as  a  discontinuous  change  in  elastic  properties  of  adja¬ 
cent  materials.  Some  of  these  aspects  shall  be  discussed  in  Sec¬ 
tion  1.4  as  they  complicate  the  models  upon  which  the  signal  processor 
is  based. 


I.  3  Basic  Wave  Characteristics 

In  the  context  of  the  idealized  assumptions  stated  in  Section 
1.1,  characteristics  of  the  basic  types  of  waves  shall  be  stated  for 
future  reference.  In  particular,  the  following  characteristics  shall 
be  considered: 

a)  Direction  of  particle  motion:  Each  wave  travels  by 
deforming  the  materials  in  particular  ways.  The  geo¬ 
phones  detect  the  deformations  along  each  of  the  sensitive 
axes  of  the  instruments.  Consequently,  it  is  important 

to  identify  the  direction  of  motion  induced  by  a  wave  as 
a  potential  classification  method. 

b)  Speed  of  the  wave  front:  Different  waves  travel  with 
different  velocities.  It  is  useful  to  establish  ordering 
relationships  among  the  wave  types. 

c)  Dispersive  waves:  The  speed  of  a  wave  can  vary  with 
frequency.  As  a  result,  the  shape  of  the  wave  train 
changes  as  the  distance  of  travel  increases  thereby  add¬ 
ing  a  significant  complication  to  the  correlation  of  outputs 
from  different  arrays. 

d)  Energy  attenuation  with  range:  The  wave  energy  dimin¬ 
ishes  with  distances  of  travel  due  to  a  variety  of  effects. 
For  example,  the  earth  acts  as  a  low  pass  filter  whose 
cutoff  frequency  diminishes  with  range.  The  manner  and 
magnitude  of  the  attenuation  helps  to  define  the  required 
signal  processing  gain. 

In  this  section,  the  four  basic  characteristics  defined  above 
shall  be  described  for  compressional ,  shear,  Rayleigh,  and  Love 
waves.  Emphasis  is  placed  on  the  characteristics  that  assist  in 
discriminating  between  wave  types. 
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1.3.1  Compressional  Waves 


Compressional  waves  are  also  referred  to  as  P- waves,  longi¬ 
tudinal  waves,  or  dilatational  waves.  Impulsive  sources  generate, 
primarily,  compressional  waves.  Shear  waves,  as  discussed  below, 
are  produced  at  interface  boundaries. 

a)  Direction  of  particle  motion:  Particle  motion  for  a  com¬ 
pressional  wave  is  always  alon  g  the  direction  of  travel 
of  the  wave.  It  consists  of  alternating  condensations 
and  rarefactions  of  the  material. 

b)  Speed  of  wave  front:  Ideally,  the  wave  propagates 
away  from  the  source  as  an  expanding  sphere,  at  least 
through  the  initial  layer.  The  velocity  of  the  wavefront 
can  be  expressed  as 

Vp  .  JhflL 

where  k  =  bulk  modulus 
U  =  shear  modulus 
c  =  density 

The  speed  of  a  compressional  wave  varies,  typically, 
from  300  to  800  m/sec.  in  sand  to  4,000  to  7,000  m/sec 
in  granite.  One  sees  that  a  variation  in  speed  of  at 
least  one  order  of  magnitude  is  possible  depending  upon 
the  type  of  material.  The  speed  also  depends  upon  the 
depth  of  the  material  and  upon  its  "lithology.''  Conse¬ 
quently,  the  speed  can  be  very  different  within  materials 
that  are  classified  as  the  same.  It  is  important  to  note 
that  this  variation  implies  that  it  is  very  unlikely  that 
propagation  speed  can  be  regarded  as  known  for  a 
deployed  surveillance  system  unless  circumstances  permit 
a  thorough  calibration. 


c)  Dispersive :  Compressional  waves  are  generally  regarded 
as  nondispersive.  Thus,  the  speed  of  the  wave  front 
does  not  change  significantly  with  frequency. 

d)  Energy  attenuation  with  ranee:  As  a  spherical  wave 
propagates  from  a  source,  the  energy  of  the  wave  is 
distributed  over  the  area  of  the  sphere.  As  the  radius 
increases,  the  energy /unit  area  must  diminish  as  the 
square  of  the  radius.  Thus,  the  amplitude  of  the  wave 
(i.e.  .  the  square-root  of  the  energy)  is  inversely  pro¬ 
portional  to  its  distance  from  the  source. 

In  addition  to  the  attentuation  of  the  wave  energy  due 
to  geometrical  spreading,  losses  also  result  from  absorp¬ 
tion  of  energy  by  the  transmission  medium.  Absorption 
losses  are  described  approximately  by  an  exponential 
function 


A  =  Aq  exp(-ar) 

where  A  represents  the  amplitude  for  a  given  wavelength, 
a  is  the  absorption  coefficient,  and  r  is  the  distance  to 
the  source.  The  absorption  coefficient  is  roughly  pro¬ 
portional  to  frequency.  Thus,  the  absorption  is  greater 
for  higher  frequencies  and  results  in  a  low-pass  character¬ 
istic  for  the  earth.  A  reasonable  value  might  be  0.25 
dB  /wavelength. 

To  indicate  the  relative  effects  of  spreading  and  absorp¬ 
tion,  Telford,  et  al.  [4]  provide  a  table  (p.  242)  that 
indicates  energy  losses  as  a  function  of  distance  and 
frequency  for  a  specific  example.  This  table  is  repro¬ 
duced  on  the  following  page  (Table  1.3-1). 


The  attenuation  with  frequency  due  to  absorption  is  severe. 
Also,  the  absorption  loss  as  distance  increases  dominates 
the  spreading  loss  except  at  very  low  frequencies.  It 
appears  that  the  energy  from  a  compressional  wave  will 
be  extremely  small  at  the  ranges  anticipated  for  the  seismic 
surveillance  system  (i.e.,  10  to  20  km). 

There  are  also  energy  losses  caused  by  partitioning  of 
the  wave  front  at  interface  boundaries.  The  description 
of  the  partitioning  depends  upon  the  elastic  properties 
of  the  adjoining  materials  and  upon  the  angle  of  incidence 
of  the  wave  upon  the  interface.  As  indicated  in  Refer¬ 
ence  [4]  (e.g..  Section  4.2.2m),  the  description  of  this 
energy  loss  is  very  complicated  but  substantial  losses 
can  result. 


♦ 


> 


1.3.2  Shear  Waves 

The  second  major  type  of  body  wave  is  the  shear  wave.  These 
waves  are  also  referred  to  as  transverse,  rotational,  or  S-waves. 

a)  Direction  of  particle  motion;  The  particle  motions  for  a 
shear  wave  are  orthogonal  to  the  direction  of  travel  of 
the  wave.  Whereas  particle  motion  for  a  compressional 
wave  is  one-dimensional,  there  are  two  degrees  of  free¬ 
dom  for  the  particle  motions  of  a  shear  wave.  The  motion 
can  be  polarized  to  be  contained  in  a  plane.  The  motions 
can  be  resolved  into  components  parallel  to  and  perpendi¬ 
cular  to  the  surface  of  the  earth.  The  components  are 
referred  to  as  SH  and  SV  waves. 

b)  Speed  of  the  wave  front:  The  velocity  of  a  shear  wave 
front  is  given  by 


V, 


V 
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Comparing  this  speed  with  the  speed  of  a  compressional 
wavefront  as  given  earlier,  it  follows  that 


IS  +  i 

0  3 


Clearly,  a  shear  wave  travels  more  slowly  than  a  com¬ 
pressional  wave.  For  most  rock  formations,  the  ratio  is 
bounded  as 


1 . 5  <  Vp/Vs  <  2.0  . 


c)  Dispersive :  Shear  waves  are  generally  regarded  as  non- 
dispersive . 

d)  Energy  attenuation  with  range:  The  remarks  given  for 
compressional  waves  apply  as  well  for  shear  waves. 

1.3.3  Rayleigh  Waves 

Rayleigh  waves  travel  along  the  free  surface  of  the  earth. 

In  many  situations  the  Rayleigh  wave  is  the  strongest  wave  gener¬ 
ated  bv  the  source.  Consequently,  this  wave  should  play  an  important 
role  for  seismic  surveillance. 

a)  Direction  of  particle  motion  :  Particle  motion  of  a  Rayleigh 
wave  is  always  in  the  vertical  direction.  Furthermore, 
the  motion  is  elliptical  and  retrograde  (i.e.,  counterclock¬ 
wise)  with  respect  to  the  direction  of  motion.  Thus  a 
Rayleigh  wave  will  be  sensed  by  both  a  horizontal  and 
a  vertical  geophone.  Because  of  the  retrograde  motion, 
a  phase  difference  of  90°  will  exist  between  the  two  geo¬ 
phones.  The  sign  of  the  difference  depends  upon  the 
orientation  of  the  horizontal  geophone  relative  to  the 
direction  of  travel  of  the  wave. 
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b)  Speed  of  wave  front:  The  speed  of  a  Rayleigh  wave  is 
less  than  for  either  type  of  body  wave.  It  is  approxi¬ 
mately  equal  to  nine-tenths  of  the  speed  of  a  shear 
wave. 

c)  Dispi  •  .  .  Generally,  the  Rayleigh-wave  velocity  varies 
with  frt  , l  -icy  (i.e.,  it  is  dispersive).  For  wavelengths 
that  are  small  compared  with  the  thickness  of  the  surface 
layer,  the  speed  is  approximately  nine-tenths  of  the 
shear  velocity  in  the  surface  layer.  For  long  wave¬ 
lengths,  the  speed  is  nine-tenths  of  ths  shear  velocity 

in  the  substratum.  The  velocity  for  medium  wavelengths 
is  intermediate  to  the  velocities  at  short  and  long 
wavelengths. 

d)  Energy  attenuation  with  range:  The  Rayleigh  wave  is 

the  strongest  wave  generated  by  a  compressional  source 

and  its  energy  attenuates  much  less  rapidly  than  for 

compressional  or  shear  waves.  Because  of  its  essentially 

two-dimensional  nature,  the  geometrical  attenuation  of  the 

amplitude  is  proportional  to  the  square-root  of  the  dis- 
_  1 

tance,  r  rather  than  r  as  for  P-  and  S-waves.  The 
Rayleigh  wave  experiences  absorption  losses,  particularly 
at  higher  frequencies.  Generally,  the  Rayleigh  wave 
is  observed  at  a  distance  from  the  source  as  a  low  fre¬ 
quency  motion  (e.g.,  less  than  10  Hz).  The  time  duration 
of  this  dispersive  wave  is  long  relative  to  compressional 
and  shear  waves.  Note,  also,  that  a  Rayleigh  wave  can 
suffer  no  partitioning  losses  since  it  is  a  surface  wave. 

The  Rayleigh  wave  may  be  the  principal  energy  source 
for  a  seismic  surveillance  system.  However ,  it  is  the 
most  slowly  moving  wave  so  its  arrival  should  be  pre¬ 
saged  by  earlier  arrivals  of  other  types  of  waves.  The 
characteristic  motion  of  the  Ravleieh  wave  should  permit 


r  ^ 


its  classification  by  the  signal  processor.  The  dis¬ 
persive  nature  of  the  wav:  may  hinder  the  use  of  corre¬ 
lation  techniques  for  measurements  obtained  from  differ- 
urrays . 

1.3.4  Love  Waves 

The  layering  of  the  earth  can  induce  a  wave-guide  effect 
when  a  low-speed  material  overlay's  a  higher-speed  substratum.  A 
surface  wave  called  a  Love  wave  is  generated. 

a)  Direction  of  particle  motion:  Particle  motion  for  a  Love 
wave  is  contained  entirely  in  the  horizontal  plane  and 
is  orthogonal  to  the  direction  of  travel  of  the  wave. 
Thus,  a  vertical  geophone  will  not  sense  the  passage  of 
a  Love  wave. 

b)  Speed  of  wave  front:  The  speed  of  a  Love  wave  is 
intermediate  between  the  shear  wave  velocity  at  the  sur¬ 
face  and  the  shear  wave  velocity  in  the  substrata.  Gen¬ 
erally,  the  Love  wave  will  have  a  greater  velocity  than 

a  Rayleigh  wave  so  the  Love  wave  will  arrive  at  a  geo¬ 
phone  before  the  Rayleigh  wave.  However,  the  energy 
of  the  Love  wave  is  usually  substantially  less  than  for 
a  Rayleigh  wave. 

c)  Dispersive :  Love  waves  are  dispersive  with  greater 
velocity  at  longer  wavelengths. 

d)  Energy  attenuation  with  range:  The  Love  wave  can  be 
regarded  as  an  SH  wave  (i.e.,  a  shear  wave  polarized 
to  the  horizontal  plane)  and  suffers  attenuation  in  the 
manner  described  in  Section  1.3.2. 


1.4 


Propagation  Paths 


Seismic  energy  can  be  propagated  from  a  source  to  a  sensor 
along  a  number  of  different  paths.  Three  basic  propagation  paths 
(i.e.,  direct ,  reflections .  and  refractions)  are  considered  in  this 
section.  Single  and  multiple  reflections  are  discussed  as  are  the 
refractions  that  generate  head  waves. 


1.4.1  Direct  Waves 

A  primary  propagation  path  for  seismic  energy  is  the  direct 
wave.  Direct  waves  are  transmitted  in  the  surface  layer  of  the 
earth  and  do  not  interact  with  any  interface  boundaries  between 
layers.  The  Rayleigh  and  Love  waves  described  in  Section  1.3  exist 
only  as  direct  waves.  However,  compressional  and  shear  waves 
can  also  travel  directly  from  the  source  to  the  sensor.  Among  the 
four  direct  waves,  the  compressional  wave  will  arrive  at  the  sensor 
first  with  the  shear  and  Love  waves  arriving  next.  The  Rayleigh 
wave  should  be  the  last  to  arrive.  Since  it  has  the  greatest  energy, 
it  generally  is  the  easiest  to  detect. 

The  direct  compressional  and  the  Love  waves  are  sensed 
only  by  the  horizontal  geophones  whereas  the  direct  shear,  if  it 
exists,  and  the  Rayleigh  waves  can  also  be  detected  by  the  vertical 
geophone.  With  two  horizontal  geophones,  the  direction  of  arrival  of 
a  wave  can  be  determined,  conceptually  at  least.  For  example,  a  com- 
pressional  wave  detected  on  one  horizontal  geophone  but  not  the  other 
obviously  defines  the  direction  of  arrival.  Arrival  at  a  45°  angle  is  her¬ 
alded  by  the  simultaneous  arrival  of  the  same  signal  at  each  geophone. 

For  a  direct  wave  traveling  with  velocity  V,  the  travel  time 
tp  from  source  to  sensor  is 


where  x  represents  the  distance  from  the  source  to  the  sensor. 
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The  time-of-arrival  of  the  same  wave  at  different  sensors  is  seen 
to  be  a  linear  function  of  the  distance  x.  As  is  discussed  below 


in  Section  1.4.2,  this  linear  characteristic  across  array  elements  can 
provide  a  mechanism  for  discriminating  between  direct  and  indirect 
waves . 


1.4.2  Reflected  Waves 

Consider  the  following  diagram  and  suppose  that  the  seismic 
source  is  located  at  S  and  that  geophones  are  located  at  Rj  and  R^. 


The  arrows  define  a  ray  path  for  the  wave  which  encounters  a 
layer  interface  boundary  at  Ij  and  1^.  Some  of  the  energy  is 
reflected  back  into  the  surface  layer.  Although  not  indicated,  some 
of  the  energy  is  refracted  into  the  second  layer.  The  energy  par¬ 
titioning  depends  upon  the  angle  of  incidence  and  upon  the  proper¬ 
ties  of  the  layers. 

The  angle  of  reflection  for  the  wave  equals  the  angle  of 
incidence.  If  the  wave  velocity  is  V,  the  layer  thickness  is  h, 
and  the  source-sensor  separation  distance  is  x,  the  travel  time 
tR  must  satisfy 


or 


T 


Thus,  the  travel  time  is  described  by  a  hyperbola  and  the  difference 
in  arrival  times  at  different  sensors  is  obtained  from  this  quadratic 
relationship . 


The  arrival  times  across  an  array  for  reflected  waves  exhibit 
a  markedly  different  characteristic  than  do  the  arrivals  of  direct 
waves.  Furthermore,  a  direct  wave  of  a  particular  type  must  arrive 
at  the  sensor  sooner  than  for  this  reflected  (i.e.,  they  are  equal 
only  for  h  =  0).  The  partitioning  of  energy  at  the  interface  boundary 
causes  the  reflected  wave  to  contain  less  energy  than  the  direct 
wave.  This  reduction  is  accentuated  by  the  fact  that  the  reflected 
wave  travels  a  greater  distance  in  the  same  layer  than  the  direct 
wave . 


Because  of  the  reflection  angle,  a  compressional  wave  gener¬ 
ally  will  be  detected  on  the  vertical  as  well  as  the  horizontal  geo¬ 
phones.  In  addition,  a  compressional  wave  is  partitioned  at  a 
boundary  into  both  reflected  compressional  and  shear  waves.  The 
reflected  shear  wave  will  arrive  later  than  the  compressional  wave 
and,  if  it  has  sufficient  energy,  be  detected  by  both  horizontal 
and  vertical  geophones. 

The  reflection  phenomena  can  be  repeated  at  the  boundaries 
of  lower  layers  with  the  refracted  waves  that  enter  each  subsequent 
layer.  Thus,  the  signals  arriving  at  the  sensor  may  be  the  result 
of  one  or  several  reflections.  The  difference  in  travel  times  may 
cause  an  overlapping  of  two  signals  leading  to  the  apparent  length¬ 
ening  of  a  single  signal  or  may  lead  to  the  appearance  of  two  dis¬ 
joint  signals.  The  time  separation  of  different  reflections  of  the 
same  signal  depends  upon  the  width  and  character  of  the  layers. 
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These  characteristics  generally  must  be  unknown  for  a  seismic  sur¬ 
veillance  system.  The  energy  contained  in  multiply  reflected  waves 
often  will  be  considerably  less  than  for  direct  or  singly  reflected 
waves  because  of  the  partitioning  and  travel  distances. 


1 .  1.3  Refracted  Waves 

Seismic  energy  can  be  either  reflected  or  refracted  at  the 
boundary  between  layers  having  different  elastic  properties.  Let¬ 
ting  V  and  represent  the  wave  velocity  in  two  adjacent  layers, 
the  angle  of  refraction  S9  is  determined  using  Snell's  law.  In 
particular,  it  is  known  that 


sin  9^  sinS., 


When  the  underlying  layer  permits  a  greater  wave  velocity 
(i.e.,  >  V^)  ,  the  refracted  wave  can  travel  along  the  interface 

at  the  velocity  V^.  This  occurs  when  the  angle  of  incidence 
satisfies 


sin  5 


c 

1 


V 

V 


1 

2 


thereby  causing  sinS^  =  1  (or  6^  =  90°).  The  refracted  wave  that 

c 

is  generated  at  this  critical  angle  9^  is  referred  to  as  a  head  wave 
From  Huygen's  principle,  the  head  wave  causes  energy  to  be  propa 
gated  in  the  overlying  layer  and  to  return  to  the  surface.  This 
situation  is  depicted  in  the  following  figure. 


r 


EARTH'S 


INTERFACE 

The  head  wave  is  significant  since  the  higher  velocity 
can  permit  this  energy  to  arrive  at  the  geophone  earlier  than  any 
direct  wave.  From  the  diagram,  it  should  be  apparent  that  head 
waves  are  not  possible  in  the  interval  SN .  For  seismic  surveillance, 
the  separation  distance  x  between  the  source  and  the  sensor  should 
be  sufficiently  large  that  SR^  >  SN  and  head  waves  may  be  a  signifi¬ 
cant  carrier  of  energy  for  the  geophones. 

Travel  time  t^  for  a  head  wave  traveling  a  horizontal  dis¬ 
tance  x  can  be  expressed  as 


where 


> 


t^  =  2hcos9^C/V^  . 

Q 

The  parameter  h  represents  the  thickness  of  the  layer  and  c^  ,  V^, 
have  been  defined  above.  The  travel  time  is  seen  to  be  a  linear 
function  of  the  separation  distance  x  with  slope  l/V^  and  intercept 
time  tj. 

To  illustrate  the  relation  between  arrival  times  for  different 
propagation  paths  of  the  same  wave,  consider  the  following  diagram. 
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DIRECT  WAVE 


For  separation  distances  less  than  x1 ,  no  head  waves  are  possible. 
Thus,  the  direct  wave  provides  the  initial  energy  received  by  the 
sensors  at  distances  x  <  x'  and  the  signal  is  augmented  at  a  later 
time  by  the  arrival  of  the  reflected  energy.  For  sensors  with  sepa¬ 
ration  distances  between  x'  and  x  ,  the  direct  wave  continues  to 

c 

be  the  first  arrival  with  the  head  and  reflected  waves  arriving 
subsequently.  The  head  wave  always  precedes  the  reflected  wave. 
For  sensors  beyond  the  crossover  distance  x^,  the  head  wave 
arrives  first  and  is  followed  by  the  direct  and  the  reflected  waves. 
At  the  greater  distances,  the  reflected  wave  assumes  a  linear  char¬ 
acteristic  (i.e.,  asymptotically)  and  must  also  arrive  after  the 
direct  wave. 


The  discussion  can  be  extended  to  consider  the  effects  of 
multiple  layers.  At  each  interface,  the  incident  wave  is  partitioned 
into  reflections  and  refractions.  Head  waves  are  possible  at  every 
boundary  at  which  a  lower  velocity  layer  overlays  a  higher  velocity 
layer.  The  refracted  waves  generated  at  the  first  layer  are  parti¬ 
tioned  at  the  next  layer.  Thus,  the  mechanism  car.  repeat  as  lone 
as  energy  remains  or  until  the  layering  ceases. 


1.5 


Error  Sources 


The  conclusions  described  above  are  obtained  through  the 
analysis  of  an  idealized  model  of  the  earth.  Realistically,  many  devi¬ 
ations  from  the  modeling  assumptions  are  possible  and  their  effects 
must  be  considered  in  any  signal  processing  design.  In  this  section, 
some  major  error  sources  are  reviewed  to  provide  an  indication  of 
their  significance  and  of  possible  compensation  mechanisms. 

1.5.1  The  Low  Velocity  Layer  (LVL) 

Generally,  the  surface  layer  of  the  earth  exhibits  propaga¬ 
tion  velocities  that  are  comparable  or  less  than  the  velocities  in 
water  (i.e.,  250  m/sec  to  1000  m/sec).  Within  this  low  velocity 
layer,  the  propagation  of  seismic  energy  is  characterized  by  several 
troublesome  properties.  First,  the  velocity  is  highly  variable  in 
this  layer  and  this  variability  can  affect  the  travel  times.  Also, 
the  width  of  the  layer  can  change  substantially  between  the  source 
and  the  sensor  (e.g.,  elevation  changes).  In  addition,  absorption 
of  seismic  energy  can  be  quite  high  thereby  greatly  reducing  the 
effective  propagation  range.  The  absorption  can  effectively  limit 
the  reflected  energy  to  be  restricted  to  nearly  vertical  paths. 

Since  the  layer  below  the  LVL  often  has  radically  different  elastic 
properties,  the  interface  will  act  as  a  strong  reflector.  Then, 
the  refracted  energy  is  less  than  at  other  interfaces  thereby 
reducing  the  energy  in  the  head  wave  or  in  multiple  reflections. 

The  energy  source  for  seismic  surveillance  generally  is 
located  on  the  earth's  surface  (or  flying  above  it).  Thus,  the 
effects  of  the  LVL  must  be  considered  and  often  will  be  the  pri¬ 
mary  physical  factor  limiting  the  effective  range  of  the  system. 
Whenever  possible,  it  will  be  desirable  to  comoensate  or  to  correct 
for  the  influences  of  the  LVL. 


1.5.  > 


Dipping  La ye  r s 


The  assumptions  that  the  layers  of  the  earth  are  horizontal  is 
an  idealization  often  violated  in  the  field.  Travel  time  characteristics 
for  reflections  from  a  dipping  subsurface  are  similar  to  but  different 
from  the  result  stated  in  Section  1.4.2.  In  particular,  the  travel  time 
is  still  a  hyperbolic  (i.e.,  quadratic)  function  of  the  separation 
distance  x,  but  the  axis  of  symmetry  is  no  longer  given  by  the 
time  axis.  That  is,  the  travel  time  is  different  for  geophones 
located  symmetrically  about  the  source.  This  effect  must  be  recog¬ 
nized  if  reflected  or  refracted  waves  are  used  in  a  localization 
scheme . 


1.5.  >  Diffractions 

Whenever  a  wave  encounters  an  irregularity  whose  radius  of 
curvature  is  less  than  the  wave  length,  the  energy  is  diffracted 
according  to  Huygen's  principle.  Thus,  energy  is  propagated  in 
all  directions  from  the  irregularity.  A  diffracted  wave  will  first 
reach  the  surface  at  a  point  directly  above  the  irregularity.  The 
travel  times  increase  with  the  distance  from  this  point  in  a  manner 
that  makes  it  indistinguishable  from  a  reflection. 

Consider  the  following  diagram: 
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The  travel  time  for  the  reflected  wave,  as  indicated  earlier,  is 
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For  the  diffracted  wave  with  S?  directly  above  the  irregularity, 
the  travel  time  is 
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Thus,  the  diffracted  wave  appears  the  same  as  the  reflected  wave 
but  the  apparent  depth  of  the  reflector  differs  by  a  factor  of  two. 

The  LVL,  dipping  layers,  and  diffraction  irregularities 
represent  three  major  sources  of  potential  errors  for  the  signal 
processor  but  many  others  are  possible.  Certainly,  inelasticities, 
nonhomogeneities,  and/or  anisotropies  will  introduce  changes  in  the 
propagation  characteristics  that  have  been  described.  However, 
even  the  consideration  of  these  three  error  sources  serves  to 
emphasize  the  difficulties  inherent  in  detecting  events  and  deter¬ 
mining  source  locations  from  the  highly  complex  signals  that  arrive 
at  the  sensors.  Fundamentally,  the  signal  processor  is  concerned 
with  detecting  events  characterized  by  low  signal-to-noise  ratios 
and  with  distinguishing  between  direct  waves,  reflected  waves, 
refracted  waves,  diffracted  waves,  and  multiples. 


1.5,4  Noise  Sources 

The  noise  background  against  which  a  source  signal  must 
be  detected  is  itself  complex.  One  can  regard  the  noise  as  being 
either  incoherent  or  coherent .  Coherent  noise  is  defined  to  repre¬ 
sent  energy  that  has  characteristics  that  will  be  recognizably 


similar  in  several  elements  of  the  sensor  array.  In  many  instances, 
this  energy  may  be  generated  by  a  specific  source  much  as  the 
energy  for  the  signal  source.  But  in  focusing  upon  a  particular 
source,  energy  from  other  potential  sources  becomes  interference 
and  is  a  nuisance  that  must  be  eliminated.  Noise  which  appears 
dissimilar  in  different  sensors  is  defined  as  incoherent  and  conforms 
with  the  more  standard  description  of  random  noise. 

The  effects  of  incoherent  noise  in  different  sensors  can  be 
reduced  simply  by  adding  sensor  outputs.  The  cancellation  effect 
is  proportional  to  the  square-root  of  the  number  N  of  sensors,  \!  N. 
The  attenuation  of  coherent  noise  can  be  accomplished  by  methods 
which  depend  either  upon  knowledge  of  its  characteristics  (i.e.  , 
filtering)  or  upon  the  directivity  of  the  source  (i.e.,  beam  forming). 
Suppression  of  both  coherent  and  incoherent  noise  is  a  basic  objec¬ 
tive  of  the  signal  processing  and  shall  be  discussed  later  as  the 
signal  processing  algorithms  are  defined. 


1.6 


Some  Basic  Modeling  Considerations 


Based  on  the  preceding  discussion,  it  is  useful  to  attempt  to 
formulate  modeling  concepts  that  can  serve  to  provide  direction  to 
the  development  of  appropriate  signal  processing  algorithms.  These 
general,  qualitative  remarks  are  appropriate  before  introducing  any 
detailed  mathematical  considerations. 


Suppose  the  source  is  basically  impulsive  (e.g.,  a  shell  blast 
or  an  artillery  recoil).  Then,  the  signal  introduced  into  the  earth 
should  have  a  short  duration  with  the  largest  amplitudes  appearing 
at  the  start  of  the  signal  and  decaying  rapidly  to  zero  thereafter. 
Thus,  the  frequency  content  of  the  input  signal  is  rich,  and  a  rela¬ 
tively  broadband  signal  should  be  expected.  From  the  preceding 
discussion,  it  is  unreasonable  to  expect  that  the  signal  received  at 
the  sensor  will  have  the  same  shape  or  character  as  the  original 
signal.  The  signal  distortion  results  for  a  variety  of  reasons. 

The  apparent  time  duration  of  a  received  signal  can  be  rela¬ 
tively  long  because  of  the  variety  of  propagation  modes  and  paths 
that  are  possible.  The  concatenation  of  direct,  reflected,  and 
refracted  waves  of  different  types  induces  a  total  signal  duration 
greatly  in  excess  of  the  primary  signal.  Further,  the  time  dura¬ 
tion  of  the  low-frequency  Rayleigh  wave  is  itself  lengthy.  Thus, 
the  observed  signal  can  appear  as  a  persistent  signal  that  is  not 
consistent  with  the  assumed  impulsive  character  of  the  source. 


For  the  remainder  of  the  discussion,  the  following  notation  is 

introduced.  Consider  the  signal  received  by  a  single  geophone  and 

denote  this  signal  as  s..  (t).  The  subscript  i  (i  =  1,2,3)  is  used  to 

denote  the  specific  component  of  the  three-dimensional  geophone. 

For  convention,  we  shall  assume  that  i  =  3  denotes  the  vertical  geo- 

t  h 

phone.  The  subscript  j  (j  =  1 . 2 . N)  denotes  the  j  element  in  an 

array  composed  of  N  three-dimensional  geophones.  The  symbol  s.(t) 

t  h 

shall  be  used  to  denote  the  three-dimensional  signal  at  the  j  element. 


The  primary  signal  will  be  regarded  as  the  signal  arriving  at  the 
geophone  as  a  result  of  a  single  propagation  mode  and  associated 
with  a  single  propagation  path.  For  example,  the  primary  signal 

for  a  Rayleigh  wave  could  be  denoted  as  s..  (t)  and  is  characterized 

R  th  ^  th 

by  its  arrival  time  t..  at  the  j  array  element  and  the  i  geophone. 

Every  primary  signal  is  assumed  to  be  identically  zero  prior  to 

arrival  time. 

The  signal  that  is  observed  through  the  output  of  a  specific 
sensor  reflects  the  arrival  of  several  primary  signals.  Because  of 
the  basic  assumptions  regarding  the  elasticity  of  the  earth,  the  pri¬ 
mary  signals  can  be  assumed  to  combine  linearly  (i.e.,  they  are 
additive).  Moreover,  a  primary  signal  s..  (t)  is  obtained  as  a  linear 
convolution  of  the  original  source  signal  with  the  impulse  response 
(i.e.  .  a  Green's  function)  of  the  earth.  The  impulse  response  will 
depend  upon  the  propagation  path,  propagation  mode,  and  distance 
for  the  reasons  given  above.  For  example,  the  energy  in  the  signal 
is  attenuated  due  to  spreading,  absorption,  and  partitioning  effects. 
Secondly,  the  earth  has  a  low  pass  filtering  effect.  Consequently, 
the  broadband  nature  of  the  original  signal  is  reduced  with  the  filter¬ 
ing  effect  increasing  as  the  source /sensor  separation  increases.  Even 
with  the  reduction  in  the  frequency  bandwidth  of  the  signal,  the  effec¬ 
tive  time  duration  of  the  primary  signal  remains  relatively  brief  (i.e., 
a  few  seconds).  Finally,  the  dispersive  nature  of  some  waves  can 
cause  fundamental  changes  in  the  character  of  the  wave  form. 

The  signal  received  along  each  axis  of  a  three-dimensional 
geophone  represents  the  projection  of  the  direction  of  the  particle 
motion  upon  the  geophone  axis.  The  signal  s.(t)  can  be  defined  in 
terms  of  the  direction  of  arrival  of  the  signal.  One  is  concerned, 
principally,  with  the  direction  of  arrival  as  defined  in  the  horizontal 
plane  (i.e.,  the  x-y  plane  as  defined  by  the  horizontal  geophones). 

It  is  this  angle  that  defines  the  direction  from  the  geophonc  to  the 
source  and  provides  the  basic  mechanism  for  localization.  The 


deflection  of  the  arrival  direction  below  the  horizontal  (i.e.  ,  the  signal 
measured  by  the  vertical  geophone)  provides  information  that  is  use¬ 
ful  for  classifying  the  type  of  wave  that  is  received  at  a  particular 

time  t... 
i) 

To  conclude  the  discussion,  it  is  useful  to  consider  the  basic 
information  to  be  produced  by  the  signal  processor.  In  essence,  the 
processing  must  provide  the  type  of  information  that  is  typically 
displayed  in  a  space-time  plot.  For  example,  the  analysis  of  seismic 
records  is  aided  by  the  simultaneous  plotting  of  the  outputs  of 
several  geophones  as  a  function  of  elapsed  time  and  distance  from 
the  sensor  to  the  source.  The  amplitude  of  the  output  can  be 
reflected  by  the  darkness  of  the  plot  at  each  point.  The  following 
figure  is  taken  from  Grant  and  West  [2]  and  provides  a  graphical 
illustration  of  some  of  the  ideas  that  have  been  introduced.  Note 
that  the  origin  of  the  plot  represents  the  location  of  an  impulsive 
source  (i.e.,  a  shot). 

The  linear  characteristic  of  the  direct  wave  and  other  more 
slowly  moving  surface  waves  is  seen  clearly.  The  appearance  of 
a  head  wave  with  its  linear  characteristic  and  crossover  point  can 
be  seen.  Also,  the  hyperbolic  characteristic  of  a  reflected  wave  is 
shown  in  this  figure.  In  addition,  a  continuously  refracted  wave 
is  observed. 

Consideration  of  this  figure  should  serve  to  re-emphasize  the 
complexity  of  the  seismic  signals  that  are  obtained  from  an  array  of 
geophones.  The  arrival  times  provide  the  information  of  basic 
interest  for  the  analysis  although  the  shortness  of  the  signal  appears 
to  permit  the  separation  of  different  propagation  modes  and  paths. 

The  availability  of  spatial  measurements  (i.e..  the  existence  of  arrays) 
obviously  is  mandatory  for  any  localization.  We  note  also  that  the 
data  presented  in  this  figure  have  been  subjected  to  considerable 
processing  and  do  not  represent  the  raw"  data  obtained  from 
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Figure  1.6-1.  (p.  112)  from  Grant  and  West  [2]. 
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the  sensors.  In  concept,  the  signal  processor  must  be  capable  of 
producing  the  type  of  information  implicit  in  this  diagram.  Clearly, 
the  existence  of  a  source  signal  is  apparent  and  the  location  of 
the  source  has  been  identified. 
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PART  II 


GENERAL  SIGNAL  PROCESSING  CONSIDERATIONS 


ABSTRACT 

The  definition  of  signal  processing  algorithms  must  be  motivated 
through  appropriate  assumptions  regarding  the  signal  and  the  output 
noise  processes.  Fundamental  properties  of  the  signals  have  been  dis¬ 
cussed  in  Part  I.  As  indicated  there,  the  signals  must  be  expected 
to  have  an  unknown  form  that  results  from  the  many  and  complicated 
propagation  paths  of  the  energy  produced  by  the  target  source.  Thus, 
the  signal  processing  algorithms  must  be  developed  using  minimal  assump¬ 
tions  regarding  the  character  of  the  signals  and  of  the  output  noise. 
Reasonable  modeling  assumptions  are  introduced  in  this  section.  These 
assumptions  are  used  to  develop  detection  and  estimation  algorithms  that 
are  compatible  and  consistent  with  the  environment  in  which  a  seismic 
surveillance  system  must  operate.  In  developing  these  algorithms, 
attention  is  directed  to  many  concerns  that  are  fundamentally  important 
to  achieve  the  most  rapidly  convergent  adaptive  array  processing  algo¬ 
rithms  . 


II.  1  Mathematical  Models  and  Problem  Definitions 


II. 1.1  The  Genera]  Problem 

To  provide  a  framework  for  the  discussion,  assume  that  a  rec¬ 
tangular  coordinate  system  is  defined  with  prescribed  origin  and  axes 
x-east,  y-north,  and  z-vertical  to  form  a  right-handed  system.  The 

location  of  each  sensor  is  defined  in  terms  of  this  system.  That  is, 

t  h 

the  position  r.  of  the  i  sensor  in  the  array  (i  =  1,2,...,N)  is  given 
as 


r.  =  x.  i  +  y.  i  +  z.  i  .  (2.1) 

—l  l  —x  l  — y  l  — : z 


Throughout,  it  is  assumed  that  there  are  N  geophones,  possibly  3-axis, 
located  essentially  on  the  surface  of  the  earth.  The  origin  is  defined 
to  be  located  on  the  earth's  surface  with  the  consequence  that  z.  =  0 
for  all  i.  The  components  x.,  y.  are  defined  by  a  distance  r^  and  an 
azimuth  angle  0. ,  as  discussed  below. 

For  the  array  dimensions  and  propagation  distances  that  are 
reasonable  for  the  seismic  surveillance  system,  the  signals  measured 
at  each  sensor  can  be  assumed  to  be  generated  by  a  plane  wave 
propagating  with  speed  c  across  the  array  in  a  direction  B,  where 


B  = 


B  i  +  B  i  +  B  i  ; 
x  —x  y  — y  z  — z 


B2  +  B2  +  B2 
x  y  z 


=  1. 


(2.2) 


The  vector  B_  represents  the  normal  to  the  wavefront  as  it  passes  the 
array.  The  form  of  the  signal  received  by  the  array  elements  cannot 
be  assumed  to  be  known.  In  comparable  elements  (e.g.,  the  vertical 
seismometers)  ,  it  is  reasonable  to  assume  that  the  basic  character  of 
the  signal  is  similar,  differences  being  caused  primarily  by  propagation 
delays  caused  by  geometric  effects.  The  form  of  the  signal  may  be 
different  along  the  sensitive  axes  of  a  3-axis  gcophonc  for  the  reasons 
discussed  in  Part  I  of  this  report  and  reconsidered  below. 
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To  simplify  the  notation,  consider  an  array  composed  entirely 

of  vertical  seismometers.  Because  of  the  common  orientation  of  the 

sensors,  the  signal,  denoted  by  s(t),  that  would  be  received  at  the 

origin  of  the  coordinate  system  is  also  received  at  each  seismometer 

at  a  suitably  delayed  time.  Specifically,  the  received  signal  at  the 
t  h 

i  sensor  located  at  r.  can  be  represented  as 


By  referring  to  (2.3),  essential  aspects  of  the  problem  can  be  identi¬ 
fied.  From  the  received  signal  one  wants  to  determine  the  direction 
3  and  the  wave  speed  c.  Both  are  unknown  and  only  the  sensor  loca¬ 
tion  r.  is  known.  Note  that  in  a  practical  system  there  generally  will 
be  uncertainty  regarding  the  exact  location  of  the  sensor. 

The  sensor  output  will  differ  from  the  signal  s.(t)  because  of 
noise  effects  n.(t).  The  noise  can  be  regarded  as  the  composite  of 
two  disparate  processes.  There  can  be  coherent  or  directional  noise 
n.  (t)  that  emanates  from  some  point  source  of  seismic  energy.  Because 
of  its  directional  character,  it  is  similar  to  the  signal  itself  and  can  be 
regarded  as  being  generated  by  a  plane  wave  having  some  speed  and 
direction.  Secondly,  there  will  be  incoherent  (i.e.,  omnidirectional) 
noise  n.  (t)  that  exists  due  to  ambient  conditions  or  sensor  electronics. 
This  noise  component  generally  will  be  described  as  a  stationary  random 
process  with  zero  mean.  The  sensor  noise  shall  be  assumed  to  be 
statistically  independent  between  sensors  of  the  array  unlike  the 
coherent  noise.  The  mathematical  assumptions  imposed  on  the  signal 
and  the  coherent  and  incoherent  noise  will  be  discussed  in  the  subse¬ 
quent  sections. 
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To  summarize,  the  i  sensor  output  shall  be  denoted  as  z.(t) 
and  it  is  assumed  to  be  the  sum  of  the  signal  and  noise. 

z^t)  =  s.(t)  +  n.(t) 

where 

n.  (t)  =  n^(t)  +  nj*(t) . 

Implicitly,  any  term  s^.  n^ ,  n^  can  be  absent  in  the  output  signal. 

This  realization  provides  the  basis  for  defining  the  problems  that  must 
be  addressed  in  developing  the  signal  processor. 

Detection  Problem 

The  signal  s.(t)  might  not  be  present  in  the  sensor  output  z.(t) 
during  any  time  interval  Tg  <  t  <  T^.  When  the  signal  is  impulse-like 

(broadband),  this  will  generally  be  the  case.  Thus,  a  time  interval 

(T^-Tg)  must  be  considered  in  formulating  the  broadband  detection 
problem. 

Estimation  Problem 

The  signal  is  characterized  by  a  direction  (3  and  a  speed  c  that 

is  assumed  to  be  common  to  the  N  sensor  outputs  when  the  signal  is 

present.  By  estimating  6,  the  direction  of  the  source  is  established 
within  the  accuracy  of  the  estimate.  The  estimate  of  the  wave  speed 
c  provides  a  basis  for  identifying  the  type  of  wave  that  has  passed  the 
array.  Further,  an  estimate  of  the  signal  itself  or  of  its  spectral  char¬ 
acteristics  provides  some  indication  of  the  nature  of  the  source. 
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(2.4) 


(2.5) 
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For  the  mathematical  analysis  of  the  detection  and  estimation  prob¬ 
lems,  mathematical  models  for  the  signal  and  coherent  noise  signals  must 
be  defined.  But  it  is  unreasonable  to  define  very  precise  models  for  the 
signal  because  of  the  complicated  nature  of  the  propagation  medium.  Con¬ 
sequently,  it  is  desirable  to  impose  the  least  rigid  assumptions  possible 

c 

about  the  signal  model.  Initially,  we  shall  assume  only  that  s.  and  n^  have 
Fourier  transforms  and  that  n?  is  a  stationary  random  process  with  mean 
zero  and  correlation  function  Rn(-u).  Other  structural  assumptions  are 
introduced  below  in  establishing  some  specific  results. 


Signals  Generated 


Seismic  Waves 


Consider  the  relation  between  the  signals  received  at  a  sensor  and 
the  physical  wave  front  that  passes  the  sensor.  Suppose  that  we  consider 
a  three-axis  seismometer  with  the  seismometers  mounted  orthogonally  with 
a  vertical  seismometer  oriented  along  the  z-axis  and  two  horizontal  seismom¬ 
eters  aligned  with  the  x  and  the  y-axes  of  the  rectangular  coordinate  system 
defined  earlier.  The  geometry  is  depicted  below. 
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The  vector  y  denotes  the  unit  vector  (i.e.,  ji  =  1)  defining 
the  normal  to  the  incoming  wave  front.  The  angle  $  defines  the  cone 
angle  measured  from  the  z-axis.  The  angle  $  is  limited  to  the  range 
0°  to  90°  since  waves  approach  the  sensor  through  the  earth.  A  surface 
wave  occurs  for  <j>  =  90°  and  a  completely  vertical  wave  occurs  for 
v  =  0°.  The  azimuth  angle  9  is  measured  as  positive  counter  clockwise 
from  the  x-axis  and  0°  £  0  360°. 

For  this  discussion,  suppose  that  the  sensors  are  aligned  with 
the  x,  y,  and  z  axes,  respectively.  Any  other  orientation  could  be 
considered  but  it  only  complicates  the  notation.  Each  sensor  is  omni¬ 
directional  relative  to  particle  motions  that  project  onto  the  direction 
of  the  sensitive  axis.  For  example,  particle  motions  can  be  sensed  by 
the  vertical  seismometer  regardless  of  the  azimuth  angle  9  as  long  as 
there  is  some  verticle  particle  motion  induced  by  the  passage  of  the 
wave  front. 

Let  us  now  consider  the  signals  induced  by  the  passage  of  P- waves, 
S-waves.  Rayleigh  waves,  and  Love  waves,  respectively. 

P-Waves 

A  compressional  wave  induces  particle  motions  along  the  direction 

of  travel  of  the  wave  front.  Suppose  that  a  direct  P-wave  travels  hori- 

'  ,y  between  the  source  and  the  sensor  along  the  x-axis.  Then,  the 

passage  of  the  wave  front  is  sensed  only  by  the  x-seismometer  and  not 

by  either  the  y-seismometer  or  by  the  vertical  seismometer.  Let  the 

particle  motion  be  denoted  by  the  signal  s(t).  In  the  absence  of  noise, 

T 

the  output  of  the  seismometer  triad  would  be  the  vector  signal  s  (t)  = 
(s(t),  0,  0),  For  an  arbitrary  P-wave,  the  direction  of  travel  induces 
output  signals  that  represent  the  projection  of  S  onto  the  coordinate 
axes.  In  terms  of  4*  and  0,  the  noise-free  output  signal  is  given  by 


cos  0  sin 
sin  ')  sin 


Sp ( t )  -  Sp ( t ) 


sp(n 


(2.6-P) 


From  (2.0-P)  it  is  upparcnt  that  the  outputs  of  a  three-axis 
seismometer  can  be  regarded  as  providing  a  direct  measurement  of 
the  wave  direction  B-  There  actually  will  be  an  ambiguity  in  the 
azimuth  angle  of  180°.  If  the  system  is  composed  only  of  vertical 
seismometers,  then  the  output  is  s(t)cos$.  Direct  P-waves  (i.e., 
i  =  9(P)  are  not  sensed  by  the  vertical  seismometer. 

S- Waves 

The  particle  motions  for  a  shear  wave  are  orthogonal  to  the 
direction  of  travel  of  the  wave.  Whereas  particle  motions  for  a  com- 
pressionai  wave  correspond  with  the  direction  of  travel,  there  are  two 
degrees  of  freedom  for  the  motions  induced  by  a  shear  wave.  These 
motions  can  be  polarized  to  be  contained  in  a  plane.  Then,  the 
motions  can  be  resolved  into  motions  parallel  to  and  perpendicular  to 
the  surface  of  the  earth  (i.e.,  the  x-y  plane). 

In  general,  the  motions  induced  by  the  passage  of  a  shear 
wave  can  be  regarded  as  more  complicated  than  the  P-wave  motions. 

Often,  the  motions  can  be  assumed  to  be  polarized  to  the  vertical 
plane  (or  to  the  horizontal  plane).  Then,  the  direction  of  the  particle 
motions  can  be  characterized  by  the  unit  vector  y  where  the  noise-free 
output  is 

Sg(t)  =  Sg(t)  I  (2.f>-S) 
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where  y  is  orthogonal  to  g 


T 

Y.  8  =  0. 

In  addition,  the  assumption  that  the  waves  are  polarized  to  the  vertical 
plane  ensures  that  the  normal  to  the  plane  defined  by  g  and  y  is  ortho¬ 
gonal  to  the  z-axis.  Using  vector  cross-product  notation,  we  have  the 
condition  that 

(x*s)Tiz  =  0 

where  1_  is  the  unit  vector  defining  the  z-axis. 

The  two  orthogonality  conditions  imply  that  0  and  <J>  can  be  deter¬ 
mined  from  y  except  for  a  180°  ambiguity  in  0.  Consequently ,  the  out¬ 
puts  of  the  seismometers  provide  a  measurement  of  the  wave  direction  g. 
Note  that  a  vertical  seismometer  will  be  sensitive  to  the  passage  of  any 
shear  wave  other  than  a  vertical  wave.  The  existence  of  a  vertical 
shear  wave  of  interest  for  a  seismic  surveillance  system  is  unlikely. 

Thus,  an  array  of  vertical  seismometers  is  sufficient  if  the  major  source 
of  energy  is  a  shear  wave.  As  indicated  above,  this  type  of  an  array 
is  insensitive  to  direct  P-waves.  In  many  instances,  the  latter  is  more 
likely  to  contain  sufficient  energy  to  be  detected  than  a  shear  wave. 

Rayleigh  Waves 

Particle  motion  of  a  Rayleigh  wave  is  always  contained  in  the 
vertical  plane.  Furthermore,  the  motion  is  elliptical  and  retrograde 
(i.e.,  counterclockwise)  with  respect  to  the  direction  of  motion  of  the 
wave.  Thus,  a  Rayleigh  wave  will  always  be  sensed  by  the  horizontal 
and  by  the  vertical  geophones.  Because  of  the  retrograde  motion,  a 
phase  difference  of  90°  will  exist  between  the  horizontal  and  the  vertical 
geophones.  The  sign  of  the  phase  difference  determines  the  direction 
of  travel  of  the  wave. 


Because  of  the  two  dimensional  nature  of  the  Rayleigh  wave, 
detection  of  a  Rayleigh  wave  provides  a  direct  measurement  of  the 
direction  of  the  source.  The  existence  of  a  Rayleigh  wave  can  be 
established  by  observing  a  90°  phase  difference  in  the  outputs  of  the 
horizontal  and  vertical  geophones.  The  direction  of  the  source  is 
determined  from  the  sign  of  the  phase  difference. 


Love  Waves 

Particle  motion  for  a  Love  wave  is  contained  entirely  in  the 
horizontal  plane  and,  being  a  shear  wave,  is  orthogonal  to  the  direc¬ 
tion  of  wave  motion.  Thus,  a  vertical  geophone  will  not  sense  the 
passage  of  a  Love  wave.  As  discussed  above  for  shear  waves,  detec¬ 
tion  of  a  Love  wave  implicitly  provides  a  measurement  of  the  direction 
of  the  wave  front  except  for  the  ambiguity  of  180°.  The  wave  motion 
is  characterized  by  the  unit  vector  where 

T 

4*  =  ° 

and 


Clearly,  >>  =  90°  and  0  =  0  ±  90G . 

The  noise-free  outputs  of  a  three-axis  seismometer  are  described 
essentially  by  the  time-varying  scalar  signal  s(t)  as  modulated  by  func¬ 
tions  of  the  direction  of  the  motion  of  the  wave  front.  One  could  con¬ 
sider  the  array  processing  problem  in  the  context  of  a  vector  output 
(i.e.,  the  three-dimensional  representation  of  the  outputs  of  the  three 
orthogonally  oriented  seismometers)  from  each  array  clement.  This 
formulation  complicates  the  analysis  and  shall  not  be  pursued  here. 
Instead,  we  shall  restrict  attention  to  an  array  of  scalar  elements 
which,  for  the  purpose  of  this  discussion,  are  assumed  to  be  vertical 
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seismometers.  The  discussion  can  be  extended  to  3-axis  geophoncs 
without  great  difficulty  by  incorporating  the  considerations  mentioned 
above.  Potentially,  the  amount  of  data  is  increased  by  a  factor  of 
three  for  a  3-axis  geophone  system.  As  this  is  an  important  considera¬ 
tion,  it  is  desirable  to  consider  in  detail  the  performance  gains  that 
may  be  obtained. 


i 


1 


r 

t 

!  > 


\ 
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II. 2 


Some  Array  Processing  Fundamentals 

A  planar  array  processing  problem  shall  be  considered  which  has 
the  following  general  description.  The  N  sensors  of  the  array  are 

t  h 

assumed  to  be  located  in  the  horizontal  (i.e.,  x-y)  plane  with  the  i 

T  T 

sensor  located  at  the  distance  r.  A  r.  (cos  0.,  sin  9.,  0)  from  the 

-l  —  i  l  l 

origin.  The  magnitude  of  the  distance  is  r.  and  the  angle  9.  is  meas¬ 
ured  counterclockwise  from  the  x-axis.  The  unit  vector  J5  defines  the 
direction  of  motion  of  the  incoming  wave  front  that  represents  the 
signal.  As  indicated  above,  this  vector  can  be  defined  in  terms  of  the 
cone  angle  Cp,  0  <  <p  <  90°,  and  the  azimuth  angle  >  0  <_  9  <_  360°. 

T 

6  6=  1. 


3  = 


The  speed  of  the  wavefront  is  represented  by  c. 

For  the  purposes  of  discussion,  it  is  useful  to  reference  the 

output  of  each  sensor  to  the  time  at  which  the  wave  passes  the  origin 

th 

of  the  coordinate  frame.  The  signal  that  appears  at  the  i  sensor 
shall  be  regarded  as  the  replica  of  the  signal  that  would  have  been 
sensed  at  the  origin  except  for  a  time  shift  x.  given  by 

,T 

6  r. 

x.  = - .  (2.7) 

l  c 

Thus,  the  noise- free  output  of  the  i**1  sensor  can  be  written  as 


s.(t)  =  s(t  +  t.);  i  =  1,2 . N  (2.8) 

where  s(t)  represents  the  time  history  of  the  signal  appearing  at 
the  origin. 


► 


By  introducing  the  time  shifts  t.  for  each  i,  i  =  1,2 . N,  the 

signals  obtained  from  each  sensor  can  be  time-aligned  since 

s.  ( t  -  t  i )  =  s(t),  i  =  1,  2, 

Then,  by  adding  the  time-delayed  versions  of  the  output,  the  signal 
is  enhanced 

N  N 

y  s.(t-T.)  =  ^  s(t)  =  n  s(t) . 

i=l  i=l 

This  simple  observation  forms  the  basis  for  array  processing. 

For  seismic  surveillance,  the  direction  and  speed  of  the  incoming 
wave  are  unknown.  Then,  one  can  select  a  value  for  3  and  for  c 
in  order  to  define  the  time-delay  t.,  i  =  1,2,...,N.  By  delaying  and 
summing  the  sensor  outputs,  a  "beam"  is  formed  in  the  selected  direc¬ 
tion  and  tests  can  be  defined  to  detect  the  existence  or  absence  of  a 
signal  emanating  from  the  direction  3.  This  concept  forms  the  basis, 
either  implicitly  or  explicitly,  for  any  array  processing  scheme  that  is 
used  for  seismic  surveillance.  In  the  following  paragraphs,  the  implica¬ 
tions  of  this  simple  model  are  explored  further  in  an  effort  to  motivate 
the  algorithmic  approach  stated  in  Part  III. 

t  h 

The  output  z.(t)  of  the  i  sensor  is  described  by  (2.4). 

Suppose  that  estimates  2,  c  of  the  signal  direction  and  the  wave 
speed  are  selected  and  used  to  define  delays  t.  .  Summing  the  delayed 
outputs,  one  obtains  an  estimate  of  the  signal  wave  from 

N 

S(t)  A  iSfE  Zi(t_  V 


i=l 

N  N 
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Clearly,  it'  t.  -  x.,  the  estimator  assumes  the  form 
1  1  i 

N 

S(t)  =  s(t)  +  ^^2  nj  (  t  — T.  )  . 

i=l 

The  coherent  component  ot  nTt)  is  assumed  io  emanate  from  a  different 
direction  than  s(t)  with  the  result  that  the  n.(t)  add  with  different  phases 
and  tend  to  cancel.  The  summation  of  the  N  sensor  outputs  tends  to 
reduce  the  influence  of  the  incoherent  noise.  Thus,  the  signal  s(t) 
should  be  enhanced  relative  to  the  noise.  On  the  other  hand,  when 
i.  t  n  ,  the  signals  also  tend  to  cancel  with  the  result  that  the  signal 
estimate  should  indicate  an  absence  of  signal  power  relative  to  the  noise. 


II. 2.1  Response  of  a  Beam-formed  Array 


The  preceding  discussion  has  shown  that  the  signal  is  enhanced 
when  its  direction  of  arrival  and  wave  speed  correspond  to  the  estimates 
of  6  and  c.  Consider  the  response  of  the  beam-formed  array  to  signals 
that  arrive  from  an  arbitrary  direction  3  with  a  speed  c.  To  simplify  this 
analysis,  only  directional  signals  shall  be  considered  and  the  response  of 
the  array  to  these  signals  shall  be  investigated. 

Suppose  that  the  sensor  outputs  s.(t)  are  delayed  by  x.  but  a 
wave  comes  from  a  direction  J3_with  speed  c.  Thus,  for  signal  alignment, 
the  outputs  would  have  to  be  delayed  by  x..  In  this  case 


N 

8(t)  =  ^  S.(t~X.) 

i=l 

N 

=  s(t+x.  -x.)  (2.10) 

i=l 

where 
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T. 

1 


T. 

1 


A  T 

A  -v  r. 


To  facilitate  the  analysis,  suppose  that  the  signal  has  a  Fourier 
transform  so  that 


s(t)  =  - \  J  S(~^  eiot 


where  S(ca)  represents  the  Fourier  transform  of  s(t).  For  a  time-delayed 
signal,  this  is  written  as 

i  r  j  oj  ( t  -  A  ) 

s  ( t  -  A. )  =  2^  J  S(w)  e  d  v  . 

—  CO 

Consequently,  it  follows  that 
N 

s  ( t )  =  s  ( t  -  r . ) 

i=l 


hf 


J-t 


=  ~  I  S(  .;)  A  (  j,  v)  eJ  1  duj 


(2. 11) 


where 


A(to,v)  A  ^  ] 
i=l 


N  .  T 
1  .  v  r 


(2. 12) 
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The  quantity  A(.:,y)  is  called  the  array  pattern  and  is  seen  to  depend 
upon  the  frequency  and  the  wave  characteristics  as  described  by  y. 

A  more  general  definition  of  the  array  pattern  will  be  introduced  subse¬ 
quently.  The  quantity 

k  A  — 

—  c 

is  called  the  wave  number  vector  and  has  the  units  of  inverse  distance. 
Note  also  that  the  wave  length  is  defined  as 


Examination  of  (2.11)  indicates  that  A(<.,v)  serves  as  a  filter  oper¬ 
ating  on  the  signal  spectrum  S(ia).  That  is,  s(t)  is  represented  as  the 
Fourier  transform  of  the  product  of  S ( u>)  and  A(w,y).  Consequently,  the 
influence  of  the  incoming  signal  on  the  array  response  is  established  by 
investigating  the  properties  of  the  array  pattern. 

If  the  signal  is  narrowband  with  center  frequency  w  ,  the  behavior 
of  A(.ic,v>)  as  a  function  of  y  determines  the  array  response.  Obviously, 
if  5/ c  =  S/c,  then  y  =  0  and  it  follows  that  for  any  frequency 

A  (  j,  0)  =  N. 

This  observation  reconfirms  the  signal  enhancement  property  of  the 
array  when  the  sensor  outputs  are  aligned.  Nonzero  values  of  v  indi¬ 
cate  the  frequency-dependent  response  of  the  array. 


II.  2.  2  Linear  Array  Response 

Much  of  the  published  literature  [e.g.,  see  References  6,  7  or 
8]  dealing  with  arrays  is  based  on  specific  array  configurations  whose 
geometry  is  chosen  to  simplify  the  analysis  and  to  gain  useful  insights 
into  the  array  response.  To  illustrate  the  frequency  response  determina¬ 
tion  in  a  closed  form,  it  is  useful  to  consider  a  linear  array  with  equi¬ 
distant  sensors  at  a  spacing  D.  Assume  that  the  first  sensor  is  located 
at  the  origin  of  the  coordinate  system  (i.e.  ,  r.  =  0)  and  suppose  that 


the  remaining  (N-l)  sensors  are  positioned  along  the  positive  x-axis. 
Then,  it  follows  that 


r. 

—l 


iD  i 
—x 


where  i_  .  denotes  the  unit  vector  defining  the  x-axis.  Suppose  further 
that  only  the  horizontal  wave  motion  is  considered  so  that 

,T 

a  r. 

—  —l 


iD  ,T 
—  0  1 
c  —  — X 


iD  a 

=  -  -  COS  0. 

c 


Similarly ,  it  follows  that 


iD 


T.  =  -  -  COS  0. 

1 

c 


For  this  planar  situation,  the  array  pattern  assumes  the  form 


N 

A(-«  ,v)  =  ^  e 

i=l 


_jiD  |co|J^  _  c^_e|  o 


Let 


n  A  D 


cos  0  cos  0 


so 


N 

A ( •  )  A  ° 

i=l 


jin 
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But  the  summation  can  bo  expressed  as 


A  (  •  ) 


I  -  e~^N  n  \ 
1  -  e_jn  / 


I  sin(N-i/2) 
!  sin ( -;/ 2) 


.  (N+ 1 )  n 
"J  2 


(2.  13) 


The  array  pattern  A(n)  is  periodic  with  period  n  =  2ir.  The 
maxima  occur  at  n  =  2n^,  n  s  0,  ±1,  ±2,  ....  The  maximum  that 
occurs  at  n  =  0  is  referred  to  as  the  mainlobe  and 


sin  (N  n / 2 ) 
si  n  ( /  2 ) 


The  other  maxima  for  n  t  0  define  the  grating  lobes  of  the  array  and 
are  undesirable  as  they  imply  that  signals  with  directions  that  generate 
these  values  of  i  are  enhanced  as  much  as  signals  with  direction  S  and 
speed  c.  The  array  pattern  vanishes  for  r)  =  2Tim/N,  m  =  ±1,  ±2,  ..., 
±(N-1)  to  form  nulls .  Waves  corresponding  to  these  nulls  are  attenuated 
totally  for  some  frequency. 

The  mainlobe  occurs  at  n  =  0  and  the  grating  lobes  are  located 
at  n  =  ±2rr,  ±4-,  ....  It  is  desirable  to  avoid  the  influence  of  the 
grating  lobes  since  signals  with  characteristics  that  yield  values  of  n 
corresponding  to  the  location  of  the  grating  lobes  are  enhanced.  The 
grating  lobes  can  be  eliminated  for  band  limited  signals  through  the 
appropriate  choice  of  the  sensor  spacing  D. 

Suppose  defines  the  highest  possible  frequency  of  the  signal. 
Then,  for  waves  such  that  c  =c,  the  variable  n  can  be  rewritten  as 


n  =  D(cos  0  -  cos  0)j/c 

£  D(cos  0  -  cos  0)  jg  /c 
=  2tt(cos  0  -  cos  0)D/Ag 

where 

Xg  4  2  tic  /  . 

To  eliminate  the  grating  lobes,  n  should  be  restricted  to  a  magnitude 
less  than  2 rr 

| n  |  <_  2tt. 

This  can  be  achieved  for  any  angles  0,0  by  choosing  the  sensor  spacing 
such  that 

D  <  XbI2. 

This  requirement  represents  the  spatial  version  of  the  well-known  sampling 
theorem.  Ag  represents  the  shortest  wavelength  present  in  the  incoming 
signal.  Thus,  there  should  be  at  least  two  sensors  contained  in  the 
interval  defined  by  A„ . 

It  is  useful  to  continue  the  analysis  for  some  examples.  Suppose 
that  a  broadside  beam  is  considered  (i.e.,  0  =  v/2)  for  D  /  XD  =  1/2. 

Then  it  follows  that 

-it  <_  n  £  it. 

For  an  end  fire  beam  (i.e.,  0  =  0),  one  sees  that 

0  <  n  <  2i. 


Nulls  in  the  array  pattern  occur  for  n  =  2nm/N;  m  =  ±1,  ±2, 
±(\'-l).  The  location  of  the  null  depends  upon  the  wavelength,  the 
array  spacing,  the  number  of  sensors,  and  the  direction  of  arrival  of 
the  wave.  For  a  broadside  beam  and  minimal  spacing  D  =  Ag/2,  the 
first  null  occurs  at  n  =  ±2~/N  so 

cos  =  2/N . 

Since  the  maximum  occurs  at  0  =  tt/2,  the  width  of  the  mainlobe  as 
measured  by  the  distance  between  the  first  nulls  is 

bw  =  2(tt/2  -  ep. 

This  provides  a  convenient  measure  of  the  beamwidth  of  the  array. 
Note  that  for  large  N,  this  beamwidth  is  given  approximately  by 

BW  <v  4/N . 


Thus,  the  beamwidth  is  proportional  to  the  number  of  sensors  for  the  fixed 
spacing  D  =  A^/2.  For  an  arbitrary  spacing  D,  the  first  null  occurs  at 

cos  °i  =  wb 

and  induces  the  beamwidth  for  large  N  of 


BW 


2A 

ND' 


1'hus,  the  beamwidth  is  a  function  of  the  length  of  the  array  rather  than 
the  number  of  sensors. 


For  endfire  beams,  the  first  null  occurs  at 


or 


2  TT 
N 


( 1  -  cos  Oj) 


2jtD 

X 


cos  0  j 
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and  the  resulting  beamwidth  is 


BW  =  20,. 

For  small  A/ND,  this  reduces  approximately  to 


Now,  the  beamwidth  is  inversely  proportional  to  the  square  root  of 
the  array  length. 

Another  measure  of  the  array  capability  is  provided  by  deter¬ 
mining  the  magnitude  of  the  sidelobes  relative  to  the  magnitude  of  the 
mainlobe.  It  can  be  shown  for  large  N  that  the  magnitude  |A  ^  |  of  the 
first  sidelobe  to  the  magnitude  |A^  |  of  the  mainlobe  is  approximately 

IA1  I  2 

*  *  • 

Thus,  the  sidelobe  magnitude  is  approximately  13.5  dB  less  than  the 
mainlobe.  To  first  order,  this  ratio  is  not  affected  by  the  array  length 
nor  by  the  number  of  sensors. 

The  linear  array  produces  an  array  pattern  that  depends  upon 
the  beam  direction  -i.  If  the  signal  direction  is  not  known,  it  may 
be  desirable  to  consider  array  configurations  that  exhibit  angular 
symmetry.  In  this  case,  circular  arrays  may  be  more  suitable.  The 
analysis  of  a  circular  array  is  considerably  more  complicated  than  for 
a  linear  array  (e.g.,  see  Ma  [7!  )  and  shall  not  be  presented  here. 
However,  for  any  array  configuration,  the  array  pattern  can  be  studied 
numerically  by  returning  to  the  defining  relation  (2.12).  For  a  small 
array  (i.e.,  N  <  20)  the  calculations  are  not  lengthy.  For  a  planar 
array,  the  graphical  presentation  of  the  numerical  results  is  more 
complicated.  We  shall  return  to  this  point  in  Fart  III. 


Let  us  summarize  and  generalize  the  important  points  of  the 

discussion  of  linear  arrays. 

1.  The  magnitude  of  the  array  pattern  provides  the  basic 
characterization  of  the  array  response.  Because  of  the 
polynomial  character  of  the  array  pattern,  it  is  useful  to 
investigate  the  maxima  and  minima  of  ]A(u>,v)  |. 

2.  The  maximum  of  |A(<o,v)  |  is  achieved  for  v  =  0.  However, 
maxima  of  equal  magnitude  can  be  experienced  for  other 
values  of  v.  The  influence  of  these  grating  lobes  can  be 
eliminated  through  the  sensor  spacing. 

3.  Grating  lobes  can  also  be  avoided  by  using  nonuniform 
spacings  of  array  elements.  The  uniform  spacing  of  elements 
is  often  introduced  more  for  the  convenience  of  the  analyst 
than  for  any  other  defensible  reason.  Non  uniform  spacings 
may  be  more  appropriate  when  some  prior  knowledge  of  beam 
direction  is  available. 

4.  The  resolution  of  the  array  is  measured  primarily  by  the 
width  of  the  main  lobe  and  by  the  ratio  of  the  mainlobe  and 
sidelobe  magnitudes.  For  the  example,  the  width  of  this  lobe  is 
inversely  proportional  to  the  length  of  the  array.  The  relative 
magnitudes  are  independent  of  the  number  of  sensors. 

5.  As  discussed  below,  the  array  resolution  can  be  modified  for 
an  array  of  fixed  dimension  through  the  introduction  of 
weights  for  the  sensor  outputs  or  through  the  nonuniform 
spacing  of  the  sensor  elements. 

6.  An  important  indicator  of  the  discrimination  capability  of  the 
array  is  obtained  by  comparing  the  magnitude  of  the  mainlobe 
with  the  sidelobe  magnitudes.  The  array  effectiveness  can 
be  regarded  as  improving  as  the  ratio  of  the  magnitude 
increases.  The  directivity  is  a  quality  that  is  intended  to 
serve  as  a  measure  of  effectiveness  and  is  defined  below. 


The  number  of  sensors  N  serves  to  introduce  limitations  regard¬ 
ing  the  processing  gain  that  is  possible  for  an  array.  The  array 
geometry  can  be  chosen  to  provide  maximal  gain  for  a  given  number 
of  sensors.  As  indicated  above,  linear  or  even  circular  arrays  with 
uniform  element  spacings  do  not  necessarily  provide  the  optimal  con¬ 
figuration.  However,  these  uniform  configurations  do  have  the 
advantages  of  convenience.  Even  for  uniform  arrays,  processing  gains 
are  possible  relative  to  the  simple  delay  and  sum  configurations  intro¬ 
duced  above  by  introducing  weights  w^  for  the  output  of  each  sensor 
element . 

For  seismic  surveillance,  the  number  of  sensors  in  an  array 
can  be  expected  to  be  relatively  small.  As  a  result,  attention  must 
be  directed  to  the  extraction  of  a  maximal  amount  of  information  from 
the  sensor  data.  This  requires  a  careful  definition  of  the  array 
performance  measure,  of  the  system  description,  and  of  the  class  of 
admissible  processing  algorithms.  It  is  to  these  questions  that  the 
following  discussion  is  addressed. 


II .  3  The  Array  Processing  Problem  and  Array  Directivity 

The  array  processing  problem  can  be  defined  in  the  context  of 
a  parameter  and  signal  estimation  problem  and  this  approach  is  taken 
and  developed  here.  Basically,  we  shall  consider  the  estimation  of  s(t), 

3  and  c  from  the  sensor  output  z^(t)  ,  i  =  1,  2,  N.  It  is  assumed 

that  the  estimator  is  restricted  to  the  general  form 

N 

s(t)  =  w.  z.(t-T.)  (2.14) 

i=l 

where  the  delays  x.  are  selected  to  estimate  the  direction  6  and  speed 
c  of  the  incoming  wavefront.  The  weights  w.  are  used  to  improve  the 
ability  of  the  array  to  discriminate  between  signals  from  nearly  equal 
direction  with  similar  speeds.  Imbedded  in  the  delays  t.  are  the  sensor 
locations.  The  array  configuration  can  also  be  adjusted  to  enhance  the 
performance  of  the  array. 

As  indicated  in  (2.14),  the  estimator  s(t)  is  formed  as  a  linear 
combination  of  the  time-delayed  sensor  outputs.  For  generality,  one 
could  consider  complex  weights  operating  on  the  complex  representation 
of  the  output  signal  in  (2.14).  This  shall  not  be  done  here  although  the 
extension  of  the  following  results  for  real  weights  and  real  output  signals 
is  simple  [e.g.,  see  Reference  8].  It  involves  primarily  a  modest  change 
in  notation. 

The  estimator  in  (2.14)  can  be  written  in.  vector  notation  as 

s (t)  =  wf  x(t)  (2.15) 

where 

T 

w  A  (wp  w^ . w.^) 

x(t)  A  (z^t-Tj),  ^2^t_T2^  ’  zN^t_TN^) 

A.  |  x  ^  ( t ) ,  x2(t),  ....  x^(t))  . 


Whet-over  possible,  summation  notation,  as  used  in  (2.14),  will  be  replaced 
with  more  compact  vector  notation  introduced  in  (2.15).  The  equivalence 
of  the  two  representations  is  apparent. 

The  sensor  output  z.(t)  is  defined  to  have  the  form 
z^t)  =  s.  (t)  +  n.  (t) . 

If  t  represents  the  actual  signal  delay,  then 

z.(t-T.)  =  s.(t-T.)  +  n.(t-T.) 

11  ii  ii 

=  s(t)  +  n.(t-t.) 

4  x.(t). 

The  collection  of  all  N  sensor  outputs  can  be  expressed  in  vector  form  as 
x(t)  =  s(t)  .1  +  n'(t)  (2.16) 

where 

iT  4  (l.  l . i)T 

and 

n'T(t)  4  (n1(t-T1),  n2(t-t2) . nN^t_TN^)  . 

Equation  (2.16)  defines  the  signal-aligned  (or  beam-formed)  statement 
of  the  problem.  Using  (2.15),  it  follows  that 

s(t)  =  wT  [ s ( t )  1_  +  n 1  ( t )  ] 

=  s(t) 

It  is  reasonable  to  expect  that  s(t)  should  equal  s(t)  in  the  absence 
of  noise  and  other  errors.  This  implies  that  the  weights  should  satisfy 
the  normalization  condition. 


N 

yi  w.  +  w^  n' (t) . 
i=l 


w. 

1 


1. 


(2.  17) 


E 


i=l 


T 

w 


1  = 


For  example,  uniform  weights,  w^  =  for  all  i ,  j ,  must  be 
w.  =  1  /N  . 

l 

We  shall  return  to  the  signal-aligned  model  (2.16)  during  the  subsequent 
model.  First,  we  shall  examine  basic  properties  of  the  array  as  implied 
by  (2.15). 


II.  3.1  The  General  Array  Pattern 

Suppose  that  the  sensor  outputs  are  delayed  by  t. 


x. 

l 


where  3  and  c  represent  estimates  of  the  wave  front  direction  and 
speed.  Then,  for  signals  with  direction  B  and  speed  c,  the  array 
output  is  given  by 


N 

s(t)  =  y'  w.s.(t-T.) 

l  i  i 

i  =  l 


where 


N 


E 


W.  s(t-T.  +  T.) 
1  1  1 


N 


i  =  l 


v 


♦ 
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Assuming,  as  above,  that  the  signal  has  a  Fourier  transform,  this 
can  be  rewritten  as 


s(t)  = 


J-  f 

2  it  J 


S(ijj)A  (uj,  v,  w)e^Jtd. 


(2.18) 


where 


A(ta,  v,  w)  = 


N 

E 

i=l 


.  T 

-jtav  r. 


w.  e 
1 


(2.  19) 


The  quantity  A(w,  v,  w)  is  defined  as  the  array  pattern  and  has  the 
same  form  introduced  in  (2.12)  with  the  addition  of  the  weighting 
vector  w.  It  should  be  emphasized  that  (2.19)  is  valid  for  any  weight¬ 
ing  vector  that  satisfies  (2.17). 

Considerable  information  is  provided  from  the  analysis  of  the 
array  pattern.  It  is  convenient  to  rewrite  the  exponent  in  (2.19)  in 
terms  of  the  angular  variables  defining  the  arrival  directions.  Note 
that  since  the  sensors  are  assumed  to  be  contained  in  the  x-y  plane 
the  sensor  location  can  be  expressed  as 


r .  =  r. 

-l  l 


j  cos  6. 
I  sin  9. 

V  o' 


Also,  the  arrival  direction  B  can  be  written  in  terms  of  the  cone  and 
azimuth  angles  as 


B  = 


ORINCON 


(,l 


Then,  we  see  that 


T 


B  3 


:a  v  r.  =  o  |  — -  r. 

c  c  /  -1 


=  u> 


-T  T 

I  l;  6  ^ 


M  -  T^ti  =  *i-MT£i 


where 


*i 


—  o  r . 

\  -  -l 


- 

=  k  6  r. 


X  = 


2ttc 


and 


k  4  2 


=  wave  number 


But  it  follows  that 
T 

3  r.  =  r.  sin  W  cos (  0- 0. )  . 

-  -l  i  i 

It  is  useful  to  consider  the  array  pattern  as  a  function  of 
the  arrival  angles  9,  ip  for  fixed  frequencies  w,  sensor  locations  r. , 
and  array  parameters  B,  c,  w.  Thus,  it  becomes  convenient  to 
define 

AN,  v,»)  A  AN,  ip) 


and  note  from  the  preceding  discussion  that 


-j  [  4'-  -  k  r.  sin<Pcos  (  A-  9. )  ] 

A  (6,</?)  =  ^  ^  w.e  11  *  (2.20) 

i=l 

A  specific  array  configuration  defines  r.  and  0..  The  weights  can  be 
determined  from  a  variety  of  arguments  such  as  presented  below. 

The  angle  'b.  is  defined  in  terms  of  r. ,  9.,  0,  (2>  and  k.  The  wave 
numbers,  k  and  k,  depend  upon  the  signal  frequency.  For  a  narrow- 
band  frequency,  the  center  frequency  represents  the  single  frequency 
of  interest.  By  assigning  values  to  all  of  these  physical  parameters, 
the  magnitude 

j  A  (  9 ,  co)  | 2  =  A(0,  <p)A*(9,  ip) 

can  be  determined  to  establish  the  array  response  to  signals  emanating 
from  directions  (0 ,  <p) .  This  computation  should  be  regarded  as  funda¬ 
mental  to  the  analysis  of  any  array  processor. 

II. 3. 2  Optimal  Array  Directivity 

Suppose  that  one  wants  the  array  to  be  as  sensitive  as  possible 
in  the  mainbeam  direction  relative  to  the  average  of  the  sensitivity  in 
all  other  directions.  This  represents  a  reasonable  objective  when  there 
is  no  information  available  regarding  the  direction  of  potential  inter¬ 
ference  sources.  One  can  formulate  this  problem  quantitatively  in  the 
following  manner. 

2 

The  array  pattern  magnitude-squared  !  A  (  0 ,  </?)  |  describes  the 
response  of  the  beam-tormed  array  to  signals  coming  from  the  direction 
0.  The  average  response  is  obtained  by  integrating  with  respect  to 
all  possible  values  of  •),  (p . 


ZtT  7T /  Z 

Average  response  4  ^  /  dO  /  |A  (  8,  <^)  J  ^g(  0 ,  </>) sinc£>d</?  .  (2.21) 


In  (2.21),  we  assume  that  the  cone  angle  (P  is  equally  likely  for 
0  <  <P  £  'f / 2  and  arrivals  from  above  the  earth  (i.e.,  tt/2  <  </?  <  it)  are 
impossible.  The  weighting  g(9 ,ip)  is  included  to  indicate  sensor 
limitations.  For  example,  a  vertical  seismometer  cannot  sense  hori¬ 
zontal  particle  motions  but  is  otherwise  isotropic.  Then,  one  might 
choose  g(9 ,  cp)  =  cos  <p  to  reflect  this  sensitivity. 

The  directivity  of  the  array  for  an  arbitrary  direction  0O ,  (fi 
is  defined  as 


D  indicates  the  sensitivity  in  the  direction  (9  ,  <P  )  relative  to  the 

o  o 

average  sensitivity  in  all  directions.  A  general  problem  can  be 
stated  in  the  foil  wing  manner. 

Maximum  Directivity:  Choose  the  weights  w.  and  the  phases 

ill.  to  maximize  D  for  the  direction  (0  ,  (p  )  . 
l  o  o 

Basically,  this  problem  statement  assumes  the  use  of  complex 

weights  with  amplitudes  w.  and  phases  I//. .  For  our  purposes,  we 

have  chosen  the  phases  '[>.  by  defining  9,  <p .  Furthermore,  we  want 

to  maximize  the  directivity  for 
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In  this  case, 


l\J  ~  ^  /x 

-j(k-k)r.sin<^cos(9-6.) 
A(0,c p)  =  )>  ]  w.e  1 


=  1  if  k  -  k 


To  proceed,  define  the  vector  b  as 


b  A  [%v  1  e 


(2.23a) 


and  let  g  be 


g  = 


jkr^sin<^cos(  0-6^) 


jkr..sin(Z>cos  ( 9-9..) 
a  N  N  |  (2.23b) 


With  these  definitions, 


A(9 ,(p)  =  b  g 


|  A  (  9,  (p) 


?  T  *  T 

!  =  (bg)  (bXg) 


.  *T  *  T, 
=  b  g  g  b 


—  b*^  G  b 


where  (  )  represents  the  complex  conjugate  of  (  ) .  Since  b  does 
not  depend  upon  0  or  (f ,  it  follows  that 


2t[  tt  /  2 


°AV  =  2^  j  d°y  )A(  9,  <^)  J  2g  (  9,  g?)sint/?d(/> 


0  0 


*T 

=  b  B  b 
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v. 


where 


nl  2 


=  ?  /■■■/ 


g(9,  </>)  Ig^g”1")  sin  (/?  d(P 


Thus,  the  directivity  for  any  0  ,  <p  has  the  form 

J  o  T  o 


D  = 


b*TG  b 


where 


and 


A  *  T 

G  4i  (7  p 

O  ~  ^0-0 


So  =  s 


(9  ,  CP  ) 

o  o 


(2.24) 


The  directivity  D  is  maximized  by  choosing  b  in  (2.24)  to 
maximize  the  ratio  of  two  quadratic  forms.  It  can  be  shown  that 
the  matrix  B  is  positive-definite.  Using  results  from  matrix  theory  l  9) 
(i.e.,  the  Rayleigh  quotient),  the  ratio  is  maximized  when  b  is  the 
eigenvector  satisfying 


G  tav, 

o  -M 


=  B  — M 


where  X^  is  the  largest  eigenvalue  of  the  regular  pencil  (A-XB). 

When  G  is  formed  as  an  outer  product  as  in  (2.24),  the  maximum 
o 

eigenvalue  is 


X 


M 


B 


-1 

So 


and  the  eigenvector  corresponding  to 


M 


is 
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b. .  =  B  *g 

-M  -o 


Since  we  have  defined  the  phase  values  ip. ,  the  preceding 
results  specialize  in  the  following  manner.  Note  that 

A  (  3 ,  g  =  wT  1  . 

—  — o  —  — 

This  implies  that 
b*TG  b  = 

=  wTOw  . 

Further,  it  follows  that 

b  =  7  w 

where  7  =  diag[e  !]  . 

Consequently , 

b*TB  b  =  wT7*TB7w  AwTCw 
T 

where  C  A  7*  B7  .  The  directivity  reduces  to 
O  w 

D  =  ~T  ~~  .  (2.25) 

w  C  w 

The  matrix  C  is  positive-definite  since  7  and  B  are  positive- 
definite.  Thus  the  preceding  results  apply  and  D  has  the  maximum 
value 


b  g  g  b 


T.  .T 
w  l  1  w 
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(2.26a) 


D 


M 


and 


-M 


T  - 1 

r  c  1  i 


(2.26b) 


Thus,  the  gains  that  maximize  the  directivity  of  the  array  in  the  direc¬ 
tion  0,  O  are  given  by  (2.26b).  This  calculation  requires  the  determina¬ 
tion  of  the  matrix  C  and  its  inverse.  For  uniform  linear  or  circular 
arrays,  the  calculation  of  B  can  be  accomplished  in  a  closed  form  [10]. 
The  inversion  of  C  can  be  computationally  burdensome,  particularly  for 
large  arrays,  but  it  is  conceptually  straightforward. 

When  there  is  no  information  regarding  the  noise  field,  the  pre¬ 
ceding  discussion  provides  a  reasonable  approach  to  choosing  the  array 
weights  w.  Simple  examples  demonstrate  that  considerable  improvements 
are  possible.  For  example,  the  directivity  of  a  uniformly-spaced  linear 
array  with  isotropic  elements  spaced  at  intervals  of  0.  425X  has  a  direc¬ 
tivity  of  12.5.  Using  optimal  weights  obtained  by  solving  (2.26b),  the 
directivity  for  this  array  is  22.0.  Thus,  D  is  increased  almost  by  a 
factor  of  two. 


II. 3. 3  Null  Placement  and  Maximum  Directivity 

The  preceding  approach  can  be  extended  to  consider  the  case  in 
which  an  interference  signal  emanates  from  a  known  direction  (  Bj. ,  V3  j ) 
with  speed  c.  To  eliminate  the  influence  of  the  interference  signal,  one 
can  choose  to  locate  a  null  of  the  array  pattern  in  the  direction  ( Oj. ,  ^)  • 

This  can  be  accomplished  with  only  a  modest  change  in  the  preceding 
result.  Let  define  the  vector  associated  with  (0^,  ( p^) 

jkr.sin()?jCOs(  0j-O^)  jkr^sint/^cos  (0r°N> 

4  !<-■  . <•  )• 

To  obtain  a  null  in  this  direction,  the  weight  vector  b  is  required  to 
be  orthogonal  to  since  the  array  pattern  magnitude  vanishes. 


0. 


(2.27) 


Now.  we  want  to  maximize  D  subject  to  the  constraint  (2.27).  Note 
that  as  many  as  (N-l)  nulls  can  be  introduced  by  defining  constraints 
having  the  form  (2.27).  This  is  obvious  since  the  N  sensors  cause 
the  weight  vector  to  be  N -dimensional.  An  optimization  problem  can 
be  defined  when  there  are  no  more  than  (N-l)  constraints. 

The  presence  of  constraints  (2.27)  reduces  the  number  of 
degrees  of  freedom  available  for  the  optimization.  In  general,  the 
imposition  of  constraints  means  that  the  maximum  directivity  is  less 
than  can  be  achieved  for  the  unconstrained  problem.  It  has  the 
effect  that  the  optimization  problem  can  be  formulated  in  a  lower¬ 
dimensional  space.  Effectively,  the  dimension  of  the  matrix  that 
must  be  inverted  in  (2.26)  is  smaller.  This  is  accomplished  by 
solving  the  M  constraints  for  M  variables  in  terms  of  the  remaining 
(N-M)  variables.  Then,  these  M  variables  are  eliminated  in  the  direc¬ 
tivity  equation  to  define  an  unconstrained  maximization  problem 
involving  (N-M)  variables.  Because  the  constraints  are  linear,  the 
nature  of  the  original  problem  is  not  altered.  One  is  still  concerned 
with  maximizing  the  ratio  of  two  quadratic  forms. 

Let  us  summarize  the  preceding  discussion. 

(1)  The  weights  of  the  estimator  (2.15)  can  be  selected 

to  maximize  the  directivity  (2.22).  The  optimal  weights 
are  obtained  from  (2.26)  and  requires  the  inversion  of 
an  (NxN),  positive-definite  matrix. 

(2)  As  we  shall  see  below,  the  general  structure  of  the 
optimal  solution  as  given  by  (2.26)  repeats  for  other 
problem  formulations. 

(3)  Nulls  can  be  imposed  in  the  direction  of  as  many  as 
(N-l)  known  interference  sources.  The  imposition  of 
these  constraints  does  not  alter  the  form  of  the  optimal 
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solution  other  than  to  reduce  the  dimensionality  of 
the  optimization  problem. 

(4)  The  problem  that  has  been  formulated  requires  almost 
no  information  regarding  the  nature  of  the  signal  and 
noise.  However,  the  calculation  of  the  average  response 
(2.21)  can  be  interpreted  as  assuming  that  the  noise  is 
uniformly  distributed  for  0  <_  9  £  2tt,  0  <_  0  £  tt/2.  None¬ 
theless,  this  approach  can  be  regarded  as  utilizing  a 
minimal  number  of  assumptions  regarding  the  structure 
of  the  signal  and  noise. 
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1 1 . 4  A  Statistical  Approach  to  Optimal  Array  Processing 


The  signal-aligned  model  of  the  problem  was  given  in  (2.16). 
This  model  shall  be  used  as  the  basis  for  the  remaining  discussion  in 
which  we  shall  introduce  statistical  descriptions  for  the  noise  to  obtain 
optimal  estimation  and  detection  algorithms. 

Consider  the  noise  n'(t).  During  any  time  interval,  it  is 
unrealistic  to  expect  that  the  output  noise  will  be  known  or  pre¬ 
dictable.  A  stochastic  model  is  indicated  and  n' (t)  shall  be  regarded 
as  a  realization  of  a  random  process.  For  a  limited  time  interval, 
it  is  reasonable  to  assume  that  n'(t)  is  a  stationary  random  process 
with  zero  mean 

E[n'  (t)  ]  =  0  for  all  t  . 


The  assumption  of  zero  mean  can  be  justified  since  the  noise  basically 
is  composed  of  a  combination  of  seismic  waves  from  many  unknown 
sources  and  electronic  sensor  noise.  If  the  sensors  are  calibrated 
adequately,  the  sensor  noise  should  have  zero  mean. 

The  second-order  statistics  for  the  noise  must  also  be  specified, 
t  h 

For  the  i  sensor,  the  correlation  of  the  output  noise  is  given  by 

r**(A)  A  E  [  n.  ( t)  n.  (t- A)  ]  . 
n  —  ii 

t  h  t  h 

The  cross-correlation  of  the  noise  in  the  i  and  j  sensors  is 


r1J(A)  A  E[n.(t)  n.(t-A)  ]  . 
n  -  l  J 

Thus,  the  correlation  matrix  for  the  noise  vector  is  an  (N”N)  matrix 

with  elements  r'^  (  A) 
n 
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Rn,(A)  =  Ef  n'  ( t)  n'T  ( t-A  )  ] 

=  {  r  ^  ( A  +  T.  -  T. )  i 

(  n  j  if 

th 

where  t.  represents  the  delay  associated  with  the  i  sensor  as 

detined  above.  It  is  assumed  that  R  ,(A)  is  either  known  or  can  be 

n 

determined.  Then,  optimal  array  processing  algorithms  can  be 
defined  for  estimating  and/or  detecting  the  signal. 

In  general,  the  noise  can  be  regarded  as  arising  from  either 
the  electronics  of  the  sensor  or  from  seismic  sources.  Electronic 
sensor  noise  can  be  expected  to  be  statistically  independent  from 
sensor  to  sensor,  thereby  inducing  a  spatial  ''white  noise."  Further¬ 
more,  it  is  reasonable  to  assume  that  this  type  of  noise  is  temporally 
white  (i.e.,  broadband  relative  to  the  range  of  frequencies  important 
for  seismic  signal  processing).  In  addition  to  the  electronic  noise, 
the  sensor  output  can  be  expected  to  reflect  noise  signals  generated 
from  unknown  seismic  sources.  Background  noise  from  a  number  of 
distant  sources  may  combine  to  produce  outputs  that  appear  to  have 
no  directional  characteristics  that  induce  a  spatial  independence  among 
the  sensor  outputs.  However,  it  is  realistic  to  expect  that  some 
directional  noise  sources  will  appear  in  the  sensor  outputs  that  will 
produce  correlation  among  the  sensor  outputs.  Commonly,  these 
directional  noise  signals  represent  the  greatest  source  of  difficulty 
for  the  signal  processor.  Their  influence  must  be  eliminated  to 
enhance  the  ability  of  the  signal  processor  to  detect  signals  ema¬ 
nating  from  specific  directions. 

Two  types  of  directional  noise  can  be  considered.  In  some 
cases  the  direction  of  the  noise  is  known.  Then,  the  signal  processor 
attempts  to  place  a  null  of  the  array  pattern  in  the  known  noise 


direction  in  order  to  eliminate  its  influence.  In  many  other  cases, 
the  existence  and  location  of  a  directional  source  is  unknown.  Then, 
the  signal  processor  should  adapt  to  recognize  and  to  eliminate  the 
effects  of  a  noise  source  having  these  unknown  characteristics. 

II.  4.1  The  Linear.  Minimum  Mean-Square  Error  Estimator 

To  begin  the  development  of  array  signal  processing  algo¬ 
rithms,  suppose  that  the  noise  covariance  matrix  R  ,  is  known. 

Later,  we  shall  consider  the  more  realistic  case  in  which  R  ,  must  be 
estimated.  Suppose  that  the  signal  s(t)  is  assumed  to  be  deterministic 
but  completely  unknown.  Let  us  assume  that  s(t)  is  to  be  estimated 
as  a  linear  combination  of  the  delayed  sensor  outputs  x(t).  Thus, 

s(t)  4  wT  x(t) 

T 

where  w  4  (w ^ ,  w^,  ....  w^)  . 

Let  us  require  that  this  estimator  is  absolutely  unbiased  (11). 
That  is,  we  want  the  conditional  expectation  to  equal  the  signal. 

E[s(t)  | s ( t) ]  =  s(t)  . 

To  achieve  this  requirement,  note  that 

E  [  s  ( t)  |  s  ( t)  ]  =  E[w^  x(t)|s(t)] 

=  wTE[x(t)  |s(t)] 

=  wT[s(t)  1  ] 


since;  E(n'(t))  =  0.  But  this  reduces  to 


E[s(t)  [s(t) 


1  N 

2>, 


i=I 


;(t)  . 


Consequently,  an  unbiased  estimator  is  obtained  by  requiring  that 


N 

2>i 


1  . 


i=l 


Observe  that  uniform  weighting 


w.  =  1/N,  i  =  1,2 . N 

i 


yields  an  unbiased  estimator  and  confirms  the  earlier  discussion. 

Suppose  that  the  weights  w  are  chosen  to  minimize  the 
mean-square  error  in  the  estimation.  That  is,  choose  w  so  that 


^AMS  =  Ef  (s(t)-s(t))  )s(t)l 


(2.28) 


is  minimized.  The  Gauss-Markov  Theorem  [e.g.,  see  Reference  11] 
asserts  that  the  absolutely  unbiased,  linear,  minimum  mean-square 
estimator  of  s  is  given  by 

T  -1 
1  R  , 

SAMS^  =  ~T  -1 
l1  R  ,1 
-  n  - 


sAMS(t)  =  -AMS  -(t)' 


(2.29) 


The  variance  0g  ot  the  error  in  this  estimator  is  given  by 

A  E( (s-s) ^ |s] 

=  (irR  _1  l)"1  • 
n 
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Let  us  compare  the  optimal  weight  w^g  with  the  weight 
that  maximizes  the  array  directivity. 


-AMS  .Tr  -1  Rn’  - 

I  Rn-  I 


-M  =  C_1i- 


Both  gains  require  the  inversion  of  an  (NxN)  matrix.  The  matrix  C 
is  related  to  the  average  of  the  magnitude-squared  array  pattern. 
Consequently,  it  permits  an  interpretation  that  has  some  characteris¬ 


tics  in  common  with  the  noise  covariance  matrix  R 


However,  the 


matrices  C  and  R  ,  are  obtained  in  very  different  ways.  The  esti¬ 
mator  (2.29)  provides  a  conceptual  flexibility  that  is  very  useful. 

In  fact,  (2.29)  provides  the  basis  for  adaptive  array  processing  in 

which  the  noise  covariance  R  .  is  estimated  from  the  sensor  data. 

n' 

By  generating  estimates  of  the  noise  covariance  for  use  in  the  signal 
estimator,  the  processor  can  be  regarded  as  adapting  to  the  current 
environment.  Thus,  it  is  useful  to  examine  the  properties  of  the 
estimator  (2.29)  in  greater  detail. 

Suppose,  first,  that  R  ,  is  a  diagonal  matrix  (i.e.,  the  noise 
is  spatially  and  temporally  white) . 


Rn,  =  diag  (rnU(0)  }  . 


Then , 


diag  fl/r  (0)} 
n 


and  it  follows  that 


,  -i  _  /_J_  _L_  i  \ 

V  l  II’  22 .  NN  ] 

\  r  r  r  / 

'  n  n  n  / 
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and 


N 


1TR','  1  =  V  1/r  11  . 

n 1  -  n 


i=l 


Consequently,  the  estimator  is  given  by 


N 


sAMS(t) 


Jr  Vx.m/r11 

2  Z-u  i  n 


s  i=  1 


From  this  expression,  it  is  obvious  that  the  output  x.  is  given  the 

greatest  weight  when  it  contains  the  least  noise  (i.e.,  r  11  <  r  ^ ,  i^i), 

n  n  1 

The  accuracy  of  the  estimator  is  dictated,  primarily  by  the 
number  N  of  sensors.  For  simplicity  suppose  that  the  noise  is  iden¬ 
tically  distributed  so  that 


11  A  2 

r  A  o 

n  =  n 


for  all  i  =  1 ,  2,  . . .  ,  N 


Then , 


N 

Vi/,2 

n 

i=  1 


-,-1 


=  a  t  N  . 
n 


(2.30) 


The  error  variance  is  inversely  proportional  to  the  number  of  sensors. 
This  dependence  is  valid  even  when  the  noise  is  not  identically  distri¬ 
buted.  Only  the  constant  of  proportionality  is  changed.  In  fact, 

2 

does  not  have  to  he  diagonal  for  the  proportionality  of  0  with 
1/N  to  be  true. 
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The  observation  that  the  error  variance  is  inversely  propor¬ 
tional  to  the  number  of  sensors  reflects  the  lack  of  any  assumptions 
regarding  the  time- variation  of  the  signal  s(t).  Generally,  noise 
effects  are  reduced  by  processing  redundant  measurements.  In  the 
case  considered  thus  far,  all  redundancy  is  provided  by  combining 
the  outputs  from  the  array  sensors  (i.e.,  spatial  redundancy). 

Without  introducing  any  structure  regarding  the  time- variation  of 
signal,  further  error  reduction  is  not  possible.  We  shall  return  to 
this  aspect  later  in  the  discussion. 

II. 4. 2  The  Maximum  Likelihood  Estimator  and  Adaptive  Arrays 

The  estimator  (2.29)  can  be  obtained  from  other  points-of- view 
that  are  useful  for  considering  other  aspects  of  the  general  problem. 
In  particular,  it  is  convenient  to  assume  that  the  noise  is  jointly 
Gaussian  with  probability  density  function 

fN(n'(t))  =  (2TT)~N/2(det  Rn,)-i  exp  M  n'T  ( t)  1nl  ( t)  }  .  (2.31) 

But  this  can  be  rewritten  in  terms  of  the  unknown  signal  using 
(  2.  16).  One  obtains 

fM(n'(t))  =  (  2tt)  ^^(detR,)  2  exp  {-\  [  x(  t) -s  ( t )  1]TR  ,  ^ 

IN  —  n  —  —  n 

(x(t)-s(t)  1]  }  .  (2. 32) 

Using  (2.32),  we  can  now  determine  the  maximum  likelihood 
estimator  of  s(t).  The  likelihood  equation  [11]  for  s(t)  is 

J-  [ x( t ) — s ( t )  l]1  R  ;'[x(t)-s(t)  11  =  0  . 

os  —  -  n  —  — 

Carrying  out  the  differentiation  with  respect  to  s,  one  obtains  the 
equation 
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0. 


I  Rn.  I  x( 0  s^jL  11 
When  solved,  the  maximum  likelihood  estimator  of  s  is  seen  to  be 


T  - 1 
1  R 

sv1T  (t)  =  — x(t) 

r  r  1  i 

—  n  - 


“  SAMS(t)  * 


(2.  33) 


Thus,  under  the  assumption  that  the  noise  is  Gaussian,  the  minimum 
mean-square  and  maximum  likelihood  estimators  of  s(t)  are  identical. 

In  (2.29)  and  in  (2.33),  it  is  assumed  that  the  noise  covariance 

matrix  R  i  is  known.  For  seismic  surveillance,  it  is  unrealistic  to 
n 

assume  that  R  ,  is  known  for  all  times.  In  fact,  the  noise  covariance 
n 

matrix  must  be  estimated  from  the  sensor  outputs.  The  repeated 
estimation  of  R^,  provides  the  basis  for  adaptive  array  processing. 

The  noise  covariance  matrix  can  be  estimated  from  samples 
of  the  outputs  of  the  N  sensors.  In  the  absence  of  a  signal 

x(t)  =  n'(t)  . 


Suppose  that  the  outputs  are  samples  at  times  t^  =  kT ,  k  =  0,1 . 

M-l.  Then,  the  maximum  likelihood  estimator  of  the  noise  covariance 
matrix  is  given  by 


M-l 

R M ,  =  jj  ]£  x(iT)  xT(iT)  (2.34) 

i  =  0 

■ir-icy  of  this  estimator  depends  upon  the  number  of  inde- 
•  -  r",:>les  used  to  form  R  The  independence  of  the  samples 

;  bv  the  choice  of  the  sampling  interval  T.  As  indicated 
•  '  1  Id]),  satisfactory  estimates  are  often  obtained 


by  using  twice  as  many  independent  samples  as  there  are  sensors 
(M  -  2N)  . 


The  estimator  R  ,  must  be  computed  in  the  absence  of  any 
signal.  In  order  to  adapt  to  a  changing  noise  environment,  the 
estimator  R^,  must  be  recomputed.  Then,  R^,  must  be  inverted  in 
order  to  determine  the  weighting  vector.  To  simplify  the  calculations, 
the  inverse  R  ,  *  can  be  computed  recursively,  thereby  permitting  a 
recursive  computation  of  the  weighting  vector.  These  recursive  rela¬ 
tions  are  obtained  using  the  Snerman-Morrison  lemma  (11]  of  matrix 
algebra  which  are  stated  below. 


The  maximum  likelihood  estimator,  using  R^,,  is  given  by 

T  ~  -1 
1 A  R  1 

SML(t)  =  jT  p  - 1  r  -(t) 
t  K  ,  i 
-  n  - 


1  P 

~T 

1  P  1 


x(t) 


where 


P  A 


J2  x(iT)xT(iT) 
.  i=0 


Let  us  denote  P  explicitly  in  terms  of  the  number  of  samples  used 
in  its  formation. 


P(k+1)  4 


J^xtiT)  xT(iT) 
i=0 

( P  “  1  (  k )  +  x(kT )  xT(kT)l 


Usinu  the  Sherman-Morrison  lemma  111],  this  becomes 
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P(k+1) 


P(k) 


(2. 34) 


r 


P(k)  x(kT)  x  (kT)  P(k) 


1  +  x  (kT)  P(k)  x(kT) 


It  is  important  to  recognize  that  no  matrix  inversion  is  required 
after  P ( 1 )  is  formed. 


The  inverse  of  the  sample  covariance  matrix  is  determined, 
except  for  the  multiplication  by  the  number  of  samples,  by  evaluating 
(2.34).  The  matrix  is  caused  to  adapt  to  the  changing  noise  environ¬ 
ment  by  continually  introducing  noise  samples,  separated  sufficiently 
in  time  to  ensure  independence.  The  adaptation  is  enhanced  by 
introducing  a  fading  memory  into  the  calculation  by  redefining 
P(k+1)  as 


P(k+1)  =  ( a  (k)  P  X(k)  +  x(kT)  xT ( kT )  ]  1 


1 


The  scalar  parameter  a(k),  0  <  a  ^  I ,  is  the  fading  factor  and  causes 
the  past  measurements  to  be  given  less  weight  when  a.  <  1  than  the 
current  measurement.  It  follows  easily,  assuming  that  a  fading  factor 
is  introduced  at  every  step  that  x[(k-l)T]  is  faded  by  a(k), 
x[  ( k- 1 ) T ]  is  faded  by  a(k-l)  a(k),  and  in  general  x[(k-  2.)T]  is 
multiplied  by  a(k)  a(k-l)  .  .  .  a(k-H-l).  Thus,  older  samples 
are  given  much  less  weight  than  the  most  recent  sample.  The  fading 
memory  maximum  likelihood  estimator  is  defined  as 


T 


and 


P(k)  x(kT)xT(k'f)P(k) 

Pa(k+l)  =  P(k) - T - = - 

1  +  x  (kT)  P(k)  x(kT) 

P(k)  A  TT^T  P“(k)  . 

—  a(k) 

With  (2.35),  a  data-directed  procedure  for  determining  the  array 
weights  has  been  defined.  The  signal  is  estimated  using  (2.35)  and  the 
array  pattern  that  results  is  obtained  by  using  (2.36)  in  (2.20).  The 
array  sensitivity  can  be  established  and  signal  estimates  are  computed. 

It  is  appropriate  now  to  consider  the  problem  of  detecting  the  presence 
of  a  signal  in  the  prescribed  beam  direction. 

II. 4. 3  The  Signal  Detection  Problem 

To  provide  a  basis  for  signal  detection,  assume  as  in  (2.31)  that 
the  noise  is  multivariate  Gaussian.  Then,  at  each  sampling  time,  the 
detection  problem  can  be  stated  in  terms  of  two  hypotheses  and  H^. 

Hq:  no  signal  present 

x(t)  =  n'(t) 

Hp  signal  present 

x(  t)  =  s  (t)  1_  +  n1  (t) 

where  n'(t)  has  the  Gaussian  distribution  given  by  (2.31).  For  the 

moment,  it  is  assumed  that  the  noise  covariance  R  i  is  known.  Also, 

n 

the  signal  is  assumed  to  be  unknown  but  nonrandom. 

The  likelihood  ratio 
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r 


is  well-known  (e.g.,  see  Reference  13)  to  provide  the  optimal  detector, 
for  example,  for  the  Neyman-Pearson  criterion.  In  (2.37),  f(xiH.)  repre¬ 
sents  the  probability  density  function  for  the  output  x  under  the  hypothe¬ 
ses  fh  ,  i  =  0,  1.  From  the  definition,  the  output  x  is  seen  to  have  a 

Gaussian  distribution  under  both  hypotheses.  The  covariance  of  x  is  the 
same  for  both  cases  and  the  distributions  differ  only  because  the  mean 
values  are  different. 

f(x|HQ)  =  (  2tt)  ^  ^(det  Rn.)  1  ^ 2  exp  i  xT  (t)  rJ'  x(  t)  j  (2.38a) 

f  ( x  |  H  ^ )  =  (  2  tt)  N^"(detRn,)  1  ^ 2  exp  j- x(  t) -s  ( t)  1]  T 

•  R^[x(t)-s(t)]j..  (2.38b) 

Substituting  (2.38)  into  (2.37)  obviously  allows  the  normalizing 
constant  ( 2  tt  )  ^^(detR^t)  to  be  ignored.  To  eliminate  the  exponen¬ 
tial,  the  logarithm  can  be  taken  without  changing  the  nature  of  the 
hypothesis  test.  Then,  one  obtains  after  straightforward  manipulations, 
the  following  test. 

H1 

s(t)  j.T  Rn,'  x(t)  <  X  (2.39) 

Ho 

where  the  choice  of  X  is  based  upon  the  desired  false  alarm  and  detec¬ 
tion  probabilities.  Equation  (2.39)  indicates  that  the  output  x(t)  is 
correlated  with  the  signal  ls(t)  after  normalizing  with  the  inverse  of 

the  noise  covariance  matrix  R  }. 

n 

There  is  a  major  problem  that  arises  in  the  application  of  (2.39) 

to  the  seismic  signal  processing  problem.  Neither  the  noise  covariance 

matrix  R  t  nor  the  signal  is  known.  As  discussed  above,  the  noise 
n  n 

covariance  R^,  can  be  estimated  using  the  sample  covariance  matrix 
obtained  from  output  samples  obtained  in  the  absence  of  any  signal. 

Next,  the  signal  s(t)  can  be  estimated  using  (2.35).  The  generalized 
likelihood  ratio  test  proceeds  by  replacing  the  unknown  parameters  in 
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the  likelihood  ratio  test  by  the  best  estimates  of  the  parameters.  In 
this  case,  (2. 39)  is  rewritten  as 


H. 


s. ..  (t)  1^  R  1  x(t)  *  X 
ML  —  n  —  < 


or 


H 


1 


s^COpPxl.)  >  i 


H 


0 


But  this  further  reduces  to 


[lTPx(t)]2  Hl 

~  ~ _  > 

<  M 


T 

l1  P  1 


H 


0 


(2. 40a) 


(2.40b) 


Thus,  one  can  use  the  left-hand  side  of  either  (2.40a)  or  (2.40b)  to 
define  the  test  statistic  for  the  hypothesis  test. 

In  (2.40b),  a  quadratic  form  is  determined  as  the  test  statistic. 

2 

This  quadratic  form  could  be  recognized  as  a  x  -variable  with  a  single 

degree  of  freedom  were  P  to  be  a  deterministic  or  known  matrix  rather 

than  a  random  variable  formed  from  the  sample  covariance  matrix.  If 

2 

P  were  known,  the  hypothesis  test  could  be  evaluated  using  the  x  ~ 
distribution  to  evaluate  the  probability  of  detection  and  probability  of 
false  alarm  for  a  prescribed  threshold  value  X.  By  ignoring  the  random 
character  of  P,  approximate  probabilities  of  detection  and  false  alarm 
could  be  determined  | 1 3 1 . 
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II. 4. 4  The  Effect  of  Signal  Models 


In  the  previous  discussion,  no  assumptions  were  made  regarding 
the  character  of  the  signal.  As  a  result,  the  signal  estimator  in  (2.35) 
is  based  entirely  on  the  aligned  outputs  of  the  N  sensors.  Further 
understanding  of  the  estimation  problem  is  obtained  by  introducing 
assumptions  regarding  the  signal  and  reconsidering  the  preceding 
analysis.  As  discussed  above,  a  seismic  surveillance  system  may  be 
required  to  detect  either  narrowband  or  wide  band  signals,  depending 
upon  the  nature  of  the  target.  Both  types  of  signals  shall  be  consid¬ 
ered  below  and  related  to  the  preceding  results. 

Narrowband  Signals:  Suppose  that  the  signal  energy  is  known 
to  be  concentrated  at  some  center  frequency  wc.  It  has  been  proven 
[e.g.,  see  References  8,  12  or  1 4 1  that  for  narrowband  signals,  the 
structure  of  the  optimal  estimator  is  given  by 

s(t)  =  w*  x(t) 

where 

w  =  3  R  }  1  . 

—  n  — 

The  parameter  8  is  a  scalar  that  follows  from  the  optimality  criterion 
(e.g.,  minimum  variance,  minimum  mean-square  error,  maximum  signal- 
to-noise  ratio).  Thus,  the  results  obtained  above  are  not  changed 
significantly  when  attempting  to  detect  the  existence  of  a  narrowband 
si  gnal . 

The  relationship  between  estimators  can  be  defined  more  precisely 
and  their  examination  provides  insight  into  the  interrelations.  Suppose 
that  S  represents  the  signal  power  in  the  narrowband  signal.  The  noise 
power  is  defined  from  the  earlier  discussion  to  be 

T  - 1 

N0  '  -  V  i  • 
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Then,  the  weight  vector  that  results  in  the  largest  output  signal-to- 
noise  ratio  is  given  by 


-SNR 


R  }  1 

n  — 


(2. 41a) 


The  weight  vector  for  the  maximum  likelihood  estimator  is  given  by 


^1L 


A 


Nq  -SNR  * 


(2.41b) 


Assuming  that  the  signal  is  narrowband  with  power  S, 
error  is  minimized  with  the  weight  vector 


s  n: 


w 

— MSE 


SNo  +  No 


HVIL ' 


the  mean-square 


(2.41c) 


As  discussed  earlier,  the  structure  of  the  estimator  is  provided  essen¬ 
tially  by  forming  the  gain  (2.41a).  Estimators  for  other  optimality 
criteria  are  obtained  by  simple  scalar  operations. 


It  is  interesting  to  note  that  the  center  frequency  of  the  narrow- 
band  signal  does  not  appear  explicitly  in  the  estimator.  To  implement 
the  MMSE  estimator,  the  signal  power  S  must  be  known.  In  many  cases, 
the  maximum  likelihood  estimator  (2.41b)  appears  to  provide  the  most 
satisfying  approach. 


Broadband  Signals:  We  have  seen  that  the  estimation  error  is 

proportional  to  the  number  of  sensors.  To  reduce  the  errors  for  a  fixed 

number  of  sensors  having  a  prescribed  geometry,  one  must  consider  the 

temporal  behavior  of  the  signal.  For  a  band  limited  signal  with  maximum 

frequency  of  f  Hz.  it  is  well  known  that  the  signal  can  be  reconstructed 
m 

from  uniformly-spaced  samples  obtained  at  a  rate  of  2fm  samples /second . 
Based  on  this  realization,  the  estimation  of  the  signal  s(t)  can  be  formu¬ 
lated  in  terms  of  an  appropriate  number  of  delayed  samples  of  the  output 
of  each  sensor  in  addition  to  the  signal-aligned  samples.  Essentially, 
the  sensor  outputs  are  filtered  before  summing  and  estimating. 


Weighted  tap-delay  lines  can  be  introduced  to  accomplish  optimal 
signal  estimation  for  broadband  signals  [  15] .  The  development  of  the 
approach  reduces  to  the  previous  discussion  through  the  appropriate 
redefinition  of  the  relevant  quantities.  This  occurs  at  the  expense  of 
a  greatly  increased  dimensionality  of  the  matrices  and  vectors  that  are 
treated.  We  shall  demonstrate  the  reformulation  of  the  problem  in 
the  subsequent  paragraphs. 

Suppose  that  L  uniformly-spaced  samples  of  each  sensor  output 
are  to  be  utilized  to  estimate  the  broadband  signal  s(t).  The  sampling 
interval  is  chosen  consistent  with  the  highest  frequency  component  of 
the  signal.  The  outputs  of  each  delay  are  multiplied  by  weights  w. 
which  are  chosen  to  obtain  a  suitable  estimator  of  the  signal.  Nota- 
tionally,  let 

x[  4  Ix^t) ,  x2(t),  ....  xN  (t)  ] 
xj  A  l  Xj(t-A)  •  x2(t-A),  ....  xN  ( t- A)  1 


-L  A  [Xi[t-(L-1)  A]  ,  x2[t-(L-l)A] . Xjvj  1 1-  (L—  1 )  A]  j 

and 


Similarly,  let 


and 


ET  T 

W.  X.  4  W  X 
i=l 


(2.42) 


It  is  apparent  that  the  signal  estimator  (2.42)  has  the  same  form 
as  introduced  in  (2.15)  at  the  beginning  of  this  discussion.  The  dimension 
associated  with  the  weight  vector  w  is  now  equal  to  LN  rather  than  N.  Thus, 
the  earlier  developments  and  remarks  can  be  applied,  virtually  without  modi¬ 
fication  except  for  matters  dealing  with  dimensionality.  For  example,  the 

noise  covariance  matrix  R  ,  is ,  now,  an  (LN*LN)  matrix.  Since  it  must  be 

n 

inverted  to  determine  the  optimal  gains,  one  can  be  faced  with  a  substantial 
computational  burden.  However,  the  uniformity  of  the  delays  causes  the 
matrix  R  ,  to  be  block  Toeplitz.  Consequently,  there  exist  efficient  numeri¬ 
cal  algorithms  that  can  be  used  to  determine  R  , 

The  Toeplitz  structure  of  the  noise  covariance  matrix  is  demon¬ 
strated  easily.  Suppose  that  n^_  represents  the  noise  in  the  output 
x ,  noting  that 


^k 


£  n  i  [  t-  ( k-1)  A]  ,  n2[t-(k-l)A] . nN  [  t-(k-l)  A]  J  . 


Then,  it  is  apparent  that 


/here 


=  R 


r^  ( £  A)  =  E  |n^[  t-(k-l)  A]  n.  ( t-(k+H-l)  A]  | 
=  E  [  n.  (kA)  n.(k+a)  a]  . 


The  covariance  matrix  depends  only  on  the  delay  difference  and  it 
follows  that  it  has  the  following  structure. 

R  =  E  |  n  nT  ]  = 

n  — 


But  this  is  a  block  Toeplitz  matrix  as  was  asserted  above. 

Because  of  the  block  Toeplitz  form  of  the  noise  covariance 
matrix,  computationally-cfficient  algorithms  exist  for  finding  the 
inverse  or  for  the  solution  of  the  normal  equations  to  determine 

the  optimal  gains.  For  example,  the  Toeplitz  matrix  inversion  algo¬ 
rithm  of  Trench  [6]  could  be  used.  Alternatively,  the  multivariable 
form  of  the  Levinson  algorithm  could  be  utilized  to  determine  the  weight 
vector  [ 17) . 
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PART  III 


GENERAL  ALGORITHMS  AND  RECOMMENDATIONS 


ABSTRACT 

The  conclusions  and  major  results  of  Parts  I  and  II  are  sum¬ 
marized  in  the  following  pages.  The  theme  of  this  presentation  has 
been  to  identify  important  physical  characteristics  of  the  signals 
generated  by  a  seismic  source.  These  characteristics  have  been  used 
to  define  general  signal  processing  algorithms  that  can  provide  an 
appropriate  framework  for  the  design  and  analysis  of  specific  algo¬ 
rithms.  Because  of  the  intrinsic  complexity  of  the  operational  environ¬ 
ment  and  of  the  limitations  on  the  deployed  system,  it  is  imperative 
that  the  greatest  processing  gain  possible  is  achieved.  The  general 
approach  that  is  presented  in  Part  III  should  provide  a  performance 
baseline  for  evaluating  any  processor.  Furthermore,  it  permits  great 
flexibility  in  the  design  of  signal  processors  for  the  operational  system 
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III.  1  General _ Algorithmic  Considerations 

In  this  section,  the  lengthy  discussions  that  have  appeared 
above  are  summarized  by  identifying  the  fundamental  computations 
that  must  be  accomplished.  First,  some  general  observations  are 
required. 

1.  An  essential  aspect  of  the  detection /estimation  problem  that 

has  been  defined  appears  in  the  estimation  of  the  noise 
covariance  matrix  (or,  equivalently,  its  inverse).  The 
estimated  covariance  matrix  must  be  repeatedly  deter¬ 
mined  from  sensor  outputs.  The  reestimation  of  provides 
the  adaptive  capability  of  the  array  processor.  This  can  be 
accomplished  in  a  variety  of  ways.  For  example  [8],  the  VVidrow's 
LMS  filter  or  the  Howells-Applebaum  adaptive  processor  are 
based,  essentially,  on  the  covariance  estimation.  As  discussed 
by  Reed,  et.  al.  (12],  more  rapid  convergence  is  achieved  by 
working  directly  with  the  covariance  matrix  and  its  inverse. 

This  increased  convergence  is  achieved  at  the  cost  of  a  greater 
computational  burden.  We  propose  the  direct  approach  here 
for  two  major  reasons. 

(i)  It  is  desirable  at  this  early  stage  of  the  development 
of  a  seismic  surveillance  system  to  obtain  some  idea 
of  the  best  convergence  that  can  be  achieved.  Also, 

(ii)  It  is  likely  that  the  arrays  used  in  an  operational 
seismic  surveillance  system  may  not  use  an  excessively 
large  number  of  sensors.  Then,  the  computational 
burden  might  not  be  excessive.  Furthermore,  with  a 
limited  number  of  sensors,  rapid  convergence  becomes 
particularly  important. 


90 


ORINCON 


In  order  to  ensure  that  adapts  to  changing  noise  conditions, 
it  is  necessary  that  new  samples  of  the  sensor  outputs  be  given 


greater  emphasis  than  outdated  samples.  While  this  effect  can 
be  assured  by  defining  a  fixed  "window"  of  samples  in  computing 
R  .  finite-memory  filters  generally  are  computationally  burdensome 
and  require  considerable  amount  of  storage.  Thus,  a  recursive 
updating  procedure  that  incorporates  a  fading  memory  may  offer 
a  significant  advantage.  Although  not  discussed  here,  many 


nonrecursive  methods  can  be  used  to  obtain  R  , 

n 

FFT  and  spectral  estimation. 


including  the 


3.  The  computation  and  display  of  the  magnitude  of  the  array  pattern 

can  provide  geometric  insights  into  the  response  of  the  processor 
for  arbitrary  array  configurations  and  geometries.  As  discussed, 
the  generalized  array  pattern  incorporates  the  weighting  vector 
of  the  signal  estimator.  Although  discussed  only  in  passing,  the 
array  pattern  can  be  modified  to  incorporate  any  directional 
limitations  of  the  sensor.  For  example,  a  vertical  seismometer 
is  not  sensitive  to  a  direct  P-wave.  The  array  processor  and 
its  attendant  pattern  should  incorporate  this  information  in 
testing  for  each  type  of  seismic  wave. 


4.  Throughout  this  discussion,  restrictive  assumptions  regarding 
the  signal  structure  have  been  carefully  avoided.  In  Section 
II.  4.  4,  the  development  was  seen  to  apply  directly  to  narrow- 
band  signals.  However,  at  the  cost  of  increased  dimensionality, 
the  results  were  extended  to  broadband  signals  through  the  use 
of  tapped-delay  lines  and  some  subsequent  notational  changes. 
Reasonable  assumptions  could  be  introduced  that  would  reduce 
the  dimension  of  the  filter  development  for  broadband  signals. 

At  this  early  stage  of  development,  more  involved  models  seem 
inappropriate.  However,  separate  processors  for  narrowband 
and  for  broadband  signals  appear  to  be  desirable. 
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5.  In  the  formulation  of  the  problem,  it  is  not  difficult  conceptually 
to  place  pattern  nulls  in  the  direction  of  known  interference 
sources.  The  computational  burden  is  increased,  however,  and 
one  must  expect  the  convergence  rate  of  the  processor  to  be 
affected.  In  the  following,  we  shall  not  explicitly  describe  the 
inclusion  of  pattern  nulls,  although  this  subject  has  been  dis¬ 
cussed  in  Part  II. 

6.  There  are  many  questions  regarding  parameter  choices  and 
processor  performance  that  can  be  resolved  primarily  through 
numerical  experimentation  using  both  simulated  and  real  data. 

In  fact,  simulations  based  on  reasonable  models  of  the  sensor 
outputs  should  be  used  for  all  preliminary  testi  ig  and  for 
development  of  the  algorithms. 

7.  As  discussed  in  Part  I  and  Section  II.  1.1,  the  components  (i.e., 
the  outputs  of  the  horizontal  and  the  vertical  seismometers)  of 
the  3-axis  geophones  can  be  used  to  classify  the  type  of  seismic 
wave  that  is  detected.  By  defining  classification  procedures  that 
are  based  on  the  physical  properties  of  the  wave  propagation,  the 
range  of  speeds  and  directions  that  must  be  considered  can  be 
limited  appi-opriately  to  reduce  the  computational  burden. 

8.  The  propagation  of  seismic  energy  from  the  source  to  an  array 
of  sensors  has  a  very  complex  description.  Several  propagation 
paths  are  possible;  the  effect  of  the  earth  upon  a  specific 
propagation  path  is  itself  very  difficult  to  describe  or  predict. 
Consequently,  it  is  generally  unrealistic  to  introduce  detailed 
mathematical  models  of  the  signal  for  the  definition  of  the  pro¬ 
cessing  algorithms.  For  this  presentation,  the  primary  assump¬ 
tions  are  the  following. 


(i)  The  elements  of  a  single  array  are  spaced  sufficiently 
close  that  the  passing  wave  can  be  regarded  as  having 
a  planar  wave  front  with  a  direction  6  and  average 
speed  c. 

(ii)  The  signal  sensed  by  similar  elements  (e.g.,  the  vertical 
seismometer)  of  the  array  is  identical  except  for  a  time 
shift  that  depends  upon  the  array  geometry  and  the 
wave  front  properties. 

(iii)  Signals  are  band  limited  and  can  be  either  narrowband 
or  broadband.  Other  than  this  qualitative  description 
of  the  signal,  mathematical  models  are  not  defined  for 
the  model. 

The  processing  gain  of  the  array  for  narrowband  signals  is 
provided  primarily  by  the  number  of  sensors  in  the  array.  For  broad¬ 
band  signals,  tapped  delay  lines  can  be  used  to  incorporate  the  redun¬ 
dancy  provided  by  successive  samples  of  the  sensor  outputs.  As  dis¬ 
cussed  in  Section  II.  4,  the  Toeplitz  form  that  obtains  using  tapped 
delay  lines  permits  the  development  of  efficient  computational  algorithms 
for  inverting  R^. 
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III.  2  Structure  of  the  General  Algorithms 


A  seismic  signal  processor  must  provide  the  following  capabilities. 

1.  Detection :  The  existence  of  a  signal  must  be  detected  from 
the  sensor  outputs.  The  performance  of  the  detector  must 
be  assessed,  typically  by  determining  the  Receiver  Operating 
Curve  (ROC)  (i.e.,  the  probability  of  detection  versus  the 
probability  of  false  alarm)  . 

2.  Estimation :  The  direction  8  and  the  speed  c  of  the  detected 
wave  front  must  be  estimated.  The  quality  of  the  estimates 
in  terms  of  bias  and  error  variance  must  also  be  determined. 

In  addition,  signal  characteristics  must  be  estimated.  For 
example,  an  estimate  of  the  center  frequency  of  a  narrowband 
signal  should  be  determined.  For  a  broadband  signal,  the 
time-of-arrival  of  the  signal  may  be  a  desirable  parameter  to 
know.  Also,  the  bandwidth  of  the  signal  should  be  estimated. 

In  order  to  discriminate  between  signals  and  noise,  it  is 
imperative  that  the  noise  covariance  matrix  is  estimated.  These 
estimates  must  be  updated  continually  from  samples  that  are 
obtained  in  the  absence  of  the  signal.  As  discussed  in  Part 
II  and  summarized  below,  the  noise  covariance  matrix  (or  its 
inverse)  provides  the  basis  for  any  signal  processing  algorithm. 

3.  Classification :  From  estimates  relating  to  the  detected  signals, 
it  is  necessary  that  the  type  of  signal  source  is  determined. 

The  features  of  the  signal  have  been  discussed  in  Part  I  and 
these  characteristics  should  provide  the  basis  for  a  classifica¬ 
tion  procedure.  The  performance  of  a  classifier  will  depend 
upon  the  quality  of  the  detection /estimation  algorithms  so  will 
not  be  addressed  herein.  It  should  be  noted  that  the  source 
of  narrowband  energy  will  often  by  moving  (e.g.,  a  truck, 
tank,  or  aircraft).  The  motion  of  the  source  will  introduce 


an  additional  property  to  the  sensor  outputs  that  may  assist 
in  the  classification  problem. 


The  approach  that  has  been  discussed  in  Part  II  and  which  is 
summarized  below  is  based  upon  the  continuing  estimation  of  the  noise 
covariance  matrix.  The  estimation  must  be  accomplished  for  time-lagged 
versions  of  the  noise-only  sensor  outputs  that  are  consistent  with  the 
range  of  delays  implied  by  the  signal  directions  and  speeds  that  are  to 
be  considered.  The  covariance  matrix  is  used  in  a  processor  that  is 
formed  as  a  concatenation  of  a  beamformer  and  a  detector.  The 
beamformer  produces  estimates  of  the  signal  parameters  that  are  used 
in  a  generalized  maximum  likelihood  detector.  Source  classification  can 
be  accomplished  following  the  detection.  Let  us  now  review  the  basic 
ingredients  of  the  proposed  signal  processor. 

Array  Outputs 

z.  (t)  ,  i  =  1,  2.  . .  .  ,  N. 

The  sensor  output  records  must  be  sampled  at  a  rate  consistent 
with  highest  frequency  in  an  anticipated  signal.  The  sampling  rate 
also  must  be  chosen  to  be  compatible  with  the  delays  introduced  in  the 
beamforming  operation  and  in  the  tapped-delay  lines  useu  for  broadband 
signals.  Because  of  the  use  of  sampled-data,  beamforming  is  accomplished 
for  a  finite  number  of  beam  angles. 

Signal  Alignment  (Beamforming) 

The  signals  are  aligned  by  delaying  the  outputs  of  sensors  at 
location  tv  by  the  amount 

6  r. 

T.  =  -  — -  . 
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The  beam  direction  is  3  and  c  represents  an  estimate  of  the  wave  speed. 
For  a  fixed  3,  several  estimates  of  c  are  introduced  whose  value  depends 
upon  the  type  of  seismic  wave  that  is  sought.  The  best  estimate  of  the 
speed  c  is  based  upon  the  detector  output  (as  discussed  below) .  The 
signal-aligned  outputs  are  denoted  as 

x.(t),  i  =  1,  2 . N 

and,  implicitly,  depend  upon  3  and  c.  For  convenience,  we  shall 
denote 

xT(t)  A  [x^t),  x 2 ( t ) . xN(t)]. 

Covariance  Estimation 

It  was  shown  in  Section  II.  4  that  optimal  signal  estimators  having 
the  general  form 

s(t)  =  wT  x(t) 

are  characterized  such  that 

w  =  y  R  1  1 
—  n  — 

where 

1T  A  (1,  1 . 1). 

Different  estimation  criteria  produce  specific  values  for  y.  Consequently 
the  inverse  of  the  noise  covariance  matrix  provides  the  fundamental 
information  required  to  estimate  s(t).  Note  that  it  is  the  inverse,  not 
Rn>  that  is  required. 

For  seismic  surveillance,  it  is  unreasonable  to  expect  that  R  ^ 
is  known  a  priori.  In  fact,  it  is  not  even  realistic  to  assume  that  the 
covariance  matrix  is  constant  during  a  long  period  of  time.  Thus,  an 
adaptive  processor  will  estimate  R  and  use  the  estimate  R  to  determine 
the  weighting  vector  for  the  array  processor.  Because  it  is  desirable 


to  ensure  that  this  estimate  is  most  sensitive  to  the  current  data,  it 
seems  reasonable  to  determine  recursively.  That  is,  new  output 
data  is  incorporated  into  the  estimator  as  a  direct  update  of  the  last 
estimate.  To  reduce  the  influence  of  data  that  may  have  been  obtained 
at  a  time  more  distant  than  is  consistent  with  the  stationarity  assumption, 
a  fading  factor  a  can  be  introduced. 

» -1 

Let  P  denote  a  matrix  that  is  proportional  to  .  Then,  P 
can  be  updated  recursively  according  to 


POO  =  Pa(k)  ,  0  <  <x(k)  <  1 


where 


Pa(k+1) 


P(k)  - 


P(k)  x(kT)  xT(kT)P(k) 
l  +  xT(kT)  P(k)  x(kT) 


The  vector  x(kT)  represents  the  output  sample  that  is  used 

-  - 1 

to  obtain  the  new  estimate  of  R  (or  R  ).  Note  that  the  covariance 

n  n 

estimate  is  obtained  from  samples  for  which  it  is  assumed  that  no 
signal  is  present.  Furthermore,  samples  used  to  estimate  R^  must  be 
obtained  at  times  that  are  separated  sufficiently  that  the  samples  can 
be  regarded  as  being  statistically  independent. 


Signal  Estimation 

Using  P,  the  maximum  likelihood  estimator  of  the  signal  is  given 
by 


sML(t)  =  — ML(^+1)  ^(t) 


where 


w^lL(k+l)  4 


1T  Pa(k+1)  1 


P  (k+1)  1  . 
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The  presence  of  a  signal  is  denoted  by  hypothesis  whereas 
its  absence  is  the  hypothesis  Hq.  The  generalized  likelihood  ratio  test 
requires  the  correlation  of  the  sensor  output  with  The  test 

can  be  written  as 

H1 

SML(t)  ^  P2t(k+1)  ^(t)  <  XM 

H0 

or,  equivalently,  as 

T  2  H. 

[  l1  Pa(k+1)  x(t)l 

it  puh)  r  ;  Xm  • 

“  Ho 

When  the  LrlS  of  the  inequality  exceeds  the  threshold  A^,  the  signal 
is  said  to  be  present.  Otherwise,  the  output  is  regarded  as  noise  only. 

Note  that  the  signal  estimate  does  not  actually  have  to  be  com¬ 
puted  to  perform  the  signal  detection  test.  Again,  this  emphasizes  the 
central  role  that  must  be  played  in  the  processor  by  the  inverse  of  the 
noise  covariance  matrix  P. 


The  performance  of  the  detection  test  is  dictated  by  the  choice 
of  the  threshold  parameter  A^.  A  precise  analysis  of  the  detector  is 
difficult  since  P  '  is  a  random  quantity.  However,  estimates  of  the 
probability  of  detection  and  false  alarm,  and,  subsequently,  the  genera¬ 
tion  of  receiver  operating  curves  (ROC)  can  be  accomplished  with  the 
introduction  of  the  following  assumptions: 

1.  Suppose  that  Pa  is  not  a  random  quantity.  That  is,  treat  this 

matrix  as  the  true  and  known  value  of  R  K 

n 

2.  Assume  that  x  is  Gaussian.  Then,  the  LHS  of  the  inequality 

2 

can  be  regarded  as  a  x  "variable.  We  shall  not  pursue  the 
analysis  of  the  detector  at  this  time.  The  general  results  of 
Reference  [13|  provide  additional  details. 


Array  Pattern 


A  basic  tool  for  the  analysis  of  an  array  configuration  and 
signal  estimator  is  provided  by  examining  the  array  pattern  A(w,v,w) 
where 


Because  of  the  importance  of  this  concept  to  the  evaluation  and  analysis 
of  a  specific  array  configuration,  we  shall  discuss  it  in  considerable 
detail. 

The  analysis  is  facilitated  by  considering  the  variables 
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Then ,  one  can  write 


A(  v  V 


N 

2>. 

i-1 


;iln*ny' 


r.  cos  0. 
1  1 

r.  sin  0. 

c  l  l 


x '  ^ 


can  be  evaluated  as  a  function  of  (n  ,  n  ) 

x  y 


The  magnitude  [Alr^ 
to  obtain  insights  regarding  the  response  for  a  prescribed  array  con¬ 
figuration  {(r.  ,0.);  i  =  1,2,...,N}  and  prescribed  sensor  weights 

{w.  j  i  =  1,2 . N  }.  The  array  pattern  magnitudes  for  a  linear  and 

for  a  circular  array  are  shown  in  Figures  III  —  1  and  III—  2 ,  respectively. 


Linear  Array 

Suppose  there  are  7  (i.e.,  N  =  7)  sensors  oriented  at  a  45° 
angle  (i.e.  ,  0.  =  45°  for  all  i)  from  the  x-axis  and  located  at  the  points 
{(i.i),  i  =  1,2,...,  7}.  Consider  equally-weighted  sensors  with  unity 
weights  (i.e.,  w.  =  1)  chosen  for  convenience.  In  this  case,  we  know 
that  the  magnitude  of  the  main  lobe  equals  the  number  of  sensors. 
Furthermore,  it  is  clear  that  the  magnitude  is  periodic  with  period  2n. 

In  Figure  III—  1 ,  a  printer  plot  of  the  magnitude  of  the  array 
pattern  is  shown.  Lines  of  equal  magnitude  are  indicated  by  the  digits 
0  through  4.  The  identification  is  the  following. 


Plot  Character 

0 

1 

2 

3 

4 

Array  Magnitude 

0.1 

2 

4 

6 

7 

The  plot  shows  the  main  lobe  and  the  vestiges  of  two  grating  lobes. 
Clearly,  the  magnitude  is  less  than  2  for  most  values  of  (nx»ri  )• 

All  of  the  nulls  are  not  shown  because  of  the  scales  that  have  been 
used.  This  is  easily  remedied  by  appropriate  scaling  changes.  Note 
also  that  the  directivity  of  the  linear  array  is  suggested  by  the  linear 
character  of  the  iso-magnitude  curves.  We  shall  return  to  this  point, 
subsequently. 


Circular  Array: 

Suppose-  that  10  (i.e.,  N  =  10)  sensors  are  located  on  a  circle 
of  radius  10  (i.e.,  r^  =  10)  at  equally  spaced  angles  of  36°  (i.e., 

-  0.  =  36°).  Assume  unity  weights  (i.e.,  w,.  =  1)  for  all  sensors. 
Therefore,  we  know  that  the  magnitude  of  the  main  lobe  is  10. 

In  Figure  III—  2 ,  a  printer  plot  of  the  magnitude  of  the  array 
pattern  is  shown.  The  pattern  is  much  more  complicated  as  indicated 
by  the  iso-magnitude  lines  for  the  portion  of  the  array  pattern  that 
is  shown.  As  indicated,  there  is  a  symmetry  along  radial  lines  sep¬ 
arated  by  36°.  The  plot  characters  and  array  magnitudes  have  the 
following  correspondence. 


Plot  Character 

0 
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1 

2 

3 

4 
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Array  Magnitude 
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10 

The  plot  indicates  the  location  of  the  main  lobe  at  (0,0)  with 

secondary  lobes  scattered  about  the  ( n  ,  n  )  plane.  Nulls  can  be 

x  y 

located  from  the  plot.  The  range  of  the  plotted  variables  can  be 
extended  further  to  locate  the  grating  lobes  and  to  further  establish 
the  response  of  the  array. 

The  variables  (nx>n  )  are  related  to  the  physical  variables 
(,;,c,Sx,Sy)  of  the  problem  in  a  straightforward  manner.  It  is 
necessary  to  invert  the  relationship  given  above  in  order  to  obtain 
greater  physical  appreciation  for  the  plots  that  have  been  presented. 

It  is  reasonable  to  consider  fixed  values  of  frequency  u)  and  wave  speed 
c  for  a  beam  direction  6.  In  particular,  assume  that  c  =  c  and  suppose 
that  the  narrowband  signal  has  center  frequency  u)  .  Then,  it  is  easy 
to  establish  a  graphically  simple  relationship  between  ( nx  >  n  )  and  the 
azimuth  and  elevation  angles  (0,4>)  of  the  interference  wave. 
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Consider  constant  values  of  0  and  $,  say  6q  anu  4>q •  as  functions  of 
(  nx '  Hy )  •  Noting  that 

0 

tan  0  = 

3x 

sin  =  ±  J  B*  +  By 

it  is  easy  to  determine  for  fixed  9q  that 

ny  =  (tan  60)nx  +  y  (Bxtan  0Q  -  0y). 

Thus,  a  constant  value  of  6  induces  a  linear  relationship  between 

n  and  n  .  In  fact,  it  follows  that  all  values  of  G  have  a  common 
x  y 

point 


For  fixed  values  Vq,  one  finds  that 


Thus,  the  lines  of  constant  are  given  as  circles  in  the  ( nx>  ny) 
The  circles  are  centered  at 


and  have  radius  | ( / c )  sin<J>g|.  Since  |sin |  <_  1,  we  see  tnat 
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Note  that  the  origin  of  the  display  depends  upon  the  beam 
direction  as  defined  by  (8x.By.).  The  frequency  and  wave  speed 
c  affect  the  origin  and  the  radius  of  the  circles  of  constant  elevation 
angles  <{>q.  Except  for  the  origin,  the  radial  lines  that  define  constant 
azimuth  angles  are  unaffected  by  signal  characteristics. 

The  circle  of  largest  radius  (i.e.,  =  90°)  relates  to  a  hori¬ 

zontal  wave.  The  combination  of  the  radius  of  this  circle,  as  defined 
by  / c  (i.e.,  essentially  the  wave  length)  with  the  array  pattern 

dimensions  provides  a  graphical  means  of  determining  the  sensor 
spacings  that  are  required  to  avoid  grating  lobes. 

Note  that  the  main  lobe  is  always  located  at  the  point  nx  =  0, 
n  =  0.  The  origin  in  Figure  III—  3  is  generally  nonzero  since  the  main 
lobe  occurs  at  0q  =  0,  Thus,  the  offset  of  the  origin  is  neces¬ 

sary  to  ensure  that  the  array  is  maximally  responsive  in  the  direction 
of  the  beam. 

Let  us  illustrate  the  use  of  the  array  pattern  magnitude  plot 
and  the  transformation  from  (lx>ny)  to  ( 0 , 4>)  by  considering  the  linear 
array  whose  magnitude  is  presented  in  Figure  III- 1 .  Attention  shall 
be  restricted  to  horizontal  waves  (i.e.,  <j>  =  90°  =  .  First,  consider 

a  broadside  beam  (i.e.,  0  —  135°  or  0  =  325°).  Then,  the  center  of 
the  circle  in  Figure  III—  3  is  given  by  (  /2wc/2c)  (1,-1) .  It  is  apparent 
that  the  circle  will  be  centered  along  the  main  lobe  ridge  emanating 
from  the  origin.  The  radius  of  the  circle  is  u>  / c.  If  a  / c  =  /2n, 
then  the  array  pattern  must  include  the  grating  lobe.  This  is  seen 
by  superimposing  the  plot  of  Figure  III- 3 ,  suitable  scaled,  upon  Figure 
III—  1 .  Since  the  sensors  have  been  assumed  to  be  separated  by  a 
distance  /2/2,  we  corroborate  the  assertion  in  Section  II.  2  that  the 
spacing  of  sensors  must  be  greater  than  X / 2  to  avoid  spatial  aliasing. 


The  argument  can  be  repeated  for  any  beam  direction  and  wave 

characteristics.  In  fact,  the  analysis  of  the  aliasing  properties  of  any 

specific  array  can  be  accomplished  by  superimposing  the  content  of 

Figure  III- 3  with  suitable  scalings  on  the  plot  of  |A  ( n  >  n  )  |.  Although 

x  y 

we  have  not  pursued  it  in  the  examples,  the  array  pattern  can  include 
any  weighting  values.  Consequently,  the  array  response  can  be  examined 
in  detail  by  using  the  graphical  tools  that  have  been  introduced.  An 
interactive  graphical  display  can  be  developed  that  can  provide  the 
basis  for  a  rapid  analysis  and  evaluation  of  the  array  performance. 

Wave  and  Source  Classification 

As  discussed  in  Part  I,  a  source  can  generate  several  types 
of  seismic  waves  which  can  travel  along  a  variety  of  propagation  paths. 
The  output  of  the  signal  detector/estimator  can  be  validated  by  compar¬ 
ing  the  response  of  the  sets  of  comparable  seismometers  to  classify  the 
wave  type.  There  must  be  a  consistency  in  the  characteristic  of  the 
detected  wave  (i.e.,  bandwidths,  particle  motions,  wave  speeds,  arrival 
directions)  to  accomplish  a  realistic  classification.  Special  purpose 
algortihms  based  on  the  wave  characteristics  discussed  in  Part  I  and 
in  Section  II.  1  can  be  defined  to  accomplish  the  classification. 

A  basic  objective  of  the  seismic  surveillance  system  is  to  establish 
the  nature  of  the  source  of  the  detected  energy  and  to  locate  and, 
possibly,  track  potential  threats.  From  the  outputs  of  a  single  array, 
it  should  be  possible  to  distinguish  between  broadband  and  narrowband 
sources.  As  has  been  discussed,  the  direction  of  the  passing  wave 
front  shall  be  estimated  and  the  estimate  can  be  expected  to  provide 
some  direction-finding  capability.  However,  the  complexity  of  the 
propagation  path  may  greatly  affect  the  direction  of  travel  of  a  wave 
front  as  it  passes  from  source  to  array. 


To  accomplish  source  localization  and  tracking,  the  output  from 
a  single  array  must  be  combined  with  other  information.  The  most 
reasonable  source  of  additional  data  is  another  array  (or  several 
arrays).  We  have  not  discussed  interarray  processing  in  this  report 
as  it  was  considered  to  be  beyond  the  scope  of  the  effort.  The  pro¬ 
cessing  of  data  from  multiple  arrays  must  be  considered  as  an  integral 
part  of  the  seismic  surveillance  system. 
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III. 3  Recommendations  and  Summary 


This  report  reflects  the  basic  recommendation  of  the  effort. 

For  the  development  of  a  useful  seismic  surveillance  system,  it  is 
imperative  that  the  signal  processor  be  developed  from  a  general 
perspective  based  upon  an  understanding  of  the  physical  problem. 

To  produce  an  operational  system  that  provides  for  the  detection  of 
target  sources  in  a  timely  manner,  the  signal  processor  must  be  able 
to  extract  maximal  information  from  the  available  data.  Furthermore, 
it  must  be  sufficiently  flexible  that  nonstandard  conditions  and  unusual 
array  configurations  can  be  accommodated  without  serious  deterioration 
in  the  quality  of  the  processor  output.  Untimely  haste  in  the  imple¬ 
mentation  of  specific  algorithms  that  are  inflexible  or  ill-suited  to  the 
physical  problem  can  only  be  counterproductive. 

To  establish  the  performance  potential  of  a  seismic  surveillance 
system,  it  is  recommended  that  the  algorithms  stated  in  Section  III. 2 
and  discussed  in  Part  II  be  implemented.  The  performance  obtained 
from  this  algorithm  should  provide  a  baseline  against  which  any 
detector  (estimator /classifier)  can  be  measured.  The  generality  of 
the  algorithm  permits  the  accomplishment  of  detailed  sensitivity  analysis 
for  important  system  parameters.  The  performance  analysis  and  sensi¬ 
tivity  assessment  can  be  facilitated  by  the  graphical  display  of  the  array 
pattern  as  discussed  in  Section  III.  2. 

The  design  and  development  of  the  signal  processor,  based  upon 
the  recommended  algorithm,  should  proceed  along  two  complementary 
paths.  On  one  hand,  a  simulation  should  be  developed  of  reasonable 
sensor  outputs.  A  suitable  simulation  can  be  easily  defined  from  the 
considerations  presented  in  Part  I.  Then,  data  can  be  produced  using 
the  simulation  that  can  be  used  to  exercise  the  signal  processing  algo¬ 
rithm.  By  processing  simulated  data,  the  performance  of  the  detector/ 
estimator /classifier  can  be  judged  by  comparing  the  outputs  with  "truth." 
The  algorithms  can  be  applied  subsequently  to  data  obtained  from  con¬ 
trolled  field  experiments.  If  there  are  substantial  inconsistencies  between 
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the  simulated  and  real  data,  the  simulation  should  be  modified  to  accord 
more  closely  with  the  physical  environment.  Thus,  the  analysis  should 
permit  feedback  between  the  artificially-produced  data  (i.e.  ,  simulation) 
and  the  experience  gained  in  the  field. 

After  the  performance  of  the  general  signal  processor  has  been 
established,  attention  should  be  directed  to  the  assessment  of  the 
computational  requirements  and  capability  of  appropriate  algorithms. 

It  is  reasonable  to  expect  that  limited  versions  of  the  general  algo¬ 
rithms  may  be  required  for  implementation  in  the  operational  system. 

An  assessment  of  possible  losses  in  performance  must  be  accomplished 
during  this  portion  of  the  effort. 

To  summarize,  this  report  is  intended  to  provide  a  perspective 
on  the  physical  problem  that  should  be  used  to  guide  the  design  and 
development  of  the  signal  processor  for  a  seismic  surveillance  system. 
Important  physical  considerations  are  discussed  in  Part  I.  In  Part  II, 
array  processing  is  reviewed  and  general  adaptive  array  processing 
algorithms  are  developed.  A  unification  of  physical  considerations  and 
signal  processing  algorithm  design  is  presented  in  this  concluding  part 
of  the  report. 
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